A new interpolated pseudodifferential preconditioner for the Helmholtz equation in heterogeneous mediathanks: Funding: The work of S. Acosta and T. Khajah was partially supported by NIH award 1R15EB035359-01A1. The work of B. Palacios was partially supported by Agencia Nacional de Investigación y Desarrollo (ANID) de Chile, Grant FONDECYT Iniciación N11220772. S. Acosta would like to thank the support provided by Texas Children’s Hospital.

Sebastian Acosta Department of Pediatrics, Baylor College of Medicine and Texas Children’s Hospital, Houston, TX, USA Tahsin Khajah Department of Mechanical Engineering, The University of Texas at Tyler, Tyler, TX, USA Benjamin Palacios Department of Mathematics, Pontificia Universidad Católica de Chile, Santiago, Chile
Abstract

This paper introduces a new pseudodifferential preconditioner for the Helmholtz equation in variable media with absorption. The pseudodifferential operator is associated with the multiplicative inverse to the symbol of the Helmholtz operator. This approach is well-suited for the intermediate and high-frequency regimes. The main novel idea for the fast evaluation of the preconditioner is to interpolate its symbol, not as a function of the (high-dimensional) phase-space variables, but as a function of the wave speed itself. Since the wave speed is a real-valued function, this approach allows us to interpolate in a univariate setting even when the original problem is posed in a multidimensional physical space. As a result, the needed number of interpolation points is small, and the interpolation coefficients can be computed using the fast Fourier transform. The overall computational complexity is log-linear with respect to the degrees of freedom as inherited from the fast Fourier transform. We present some numerical experiments to illustrate the effectiveness of the preconditioner to solve the discrete Helmholtz equation using the GMRES iterative method. The implementation of an absorbing layer for scattering problems using a complex-valued wave speed is also developed. Limitations and possible extensions are also discussed.

Keywords— Wave propagation, acoustics, high frequency, pseudodifferential calculus, preconditioner, log-linear complexity

1 Introduction

For many scientific, engineering, and biomedical applications, the numerical solution of the Helmholtz equation is essential. We are particularly motivated by ultrasound-based techniques to medical therapeutics and diagnosis where computational simulations of time-harmonic wave fields are playing an increasingly important role [2, 3, 36, 44]. Unfortunately, the discretization of the Helmholtz equation by finite difference or finite element methods renders matrices that are highly indefinite, ill-conditioned and notoriously difficult to invert, especially at high frequencies [23, 19]. Mitigating these difficulties continues to be an active area of scientific research. Much of the effort has focused on designing effective and efficient preconditioners to improve the convergence of iterative solvers. Influential work includes the shifted-Laplacian preconditioner [18, 17, 10], the analytic incomplete LU preconditioner by Gander and Nataf [21, 22], the sweeping preconditioners by Engquist and Ying [16, 15] and related approaches based on domain decomposition and Schwarz-type methods [35, 37, 9, 20, 45, 38, 46, 40, 26]. See a comprehensive review by Gander and Zhang [23]. Other ingenious approaches include the use of controllability methods to obtain periodic solutions to the wave equation developed by Bristeau et al. [8, 25], Heikkola et al. [29], Grote and Tang [28], Appelo et al. [4], and a recent time-domain preconditioner by Stolk [39]. Also, there has been an effort to derive coercive or sign-definite formulations of the Helmholtz equation [34, 24].

In the present paper, we consider a Helmholtz equation of the form,

c2Δu+ω2u+iωau=f,superscript𝑐2Δ𝑢superscript𝜔2𝑢𝑖𝜔𝑎𝑢𝑓\displaystyle c^{2}\Delta u+\omega^{2}\,u+i\omega a\,u=f,italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_u + italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_u + italic_i italic_ω italic_a italic_u = italic_f , (1)

commonly employed to model ultrasound waves with angular frequency ω𝜔\omegaitalic_ω, excited by a source f𝑓fitalic_f in biological media with variable wave speed c𝑐citalic_c and damping coefficient a𝑎aitalic_a. We re-purpose the pseudodifferential analysis of the Helmholtz equation to design a semi-analytical matrix-free preconditioner in the form of a pseudodifferential operator

(Qv)(x)=(2π)d/2q(x,ξ)v^(ξ)eixξ𝑑ξ.𝑄𝑣𝑥superscript2𝜋𝑑2𝑞𝑥𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle(Qv)(x)=(2\pi)^{-d/2}\int q(x,\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xi.( italic_Q italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_q ( italic_x , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (2)

Here xd𝑥superscript𝑑x\in\mathbb{R}^{d}italic_x ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT is the spatial variable, ξd𝜉superscript𝑑\xi\in\mathbb{R}^{d}italic_ξ ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT is the frequency (Fourier dual) variable, and v^^𝑣\hat{v}over^ start_ARG italic_v end_ARG is the Fourier transform of v𝑣vitalic_v. The symbol q𝑞qitalic_q is defined in order to approximately invert the Helmholtz operator at high frequencies. The major drawback of a formulation such as (2) is that its direct evaluation is computationally expensive. Our approach to reduce the complexity of evaluating (2) is inspired by the work of Bao and Symes [7] who devised an algorithm to evaluate each term of a “classical” pseudodifferential expansion. If the symbol q𝑞qitalic_q is a term of a classical expansion, then by definition q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) is homogeneous in ξ𝜉\xiitalic_ξ [41, 27], which implies that for each fixed x𝑥xitalic_x, the function q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) is fully characterized by its profile on the unit-sphere ξ/|ξ|𝜉𝜉\xi/|\xi|italic_ξ / | italic_ξ | in the frequency domain. Consequently, Bao and Symes based their algorithm on a Fourier series approximation of q(x,ξ/|ξ|)𝑞𝑥𝜉𝜉q(x,\xi/|\xi|)italic_q ( italic_x , italic_ξ / | italic_ξ | ) which allowed them to separate the dependence on the spatial variable x𝑥xitalic_x from the frequency variable ξ𝜉\xiitalic_ξ in order to evaluate (2) efficiently. When q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) is not a classical term, then the separability between x𝑥xitalic_x and ξ𝜉\xiitalic_ξ becomes much more cumbersome to realize, especially for three-dimensional scenarios, d=3𝑑3d=3italic_d = 3. In general, it is required to sample the frequency variable ξd𝜉superscript𝑑\xi\in\mathbb{R}^{d}italic_ξ ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT to define interpolation points and their special-function interpolants to establish the desired compression [12, 13].

By contrast, we do not assume that q𝑞qitalic_q is a general symbol or a term of a classical expansion, but instead assume that

q(x,ξ)=q(c(x),ξ)𝑞𝑥𝜉𝑞𝑐𝑥𝜉\displaystyle q(x,\xi)=q(c(x),\xi)italic_q ( italic_x , italic_ξ ) = italic_q ( italic_c ( italic_x ) , italic_ξ ) (3)

for a given function c:d[cmin,cmax]:𝑐superscript𝑑subscript𝑐minsubscript𝑐maxc:\mathbb{R}^{d}\to[c_{\rm min},c_{\rm max}]\subset\mathbb{R}italic_c : blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT → [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] ⊂ blackboard_R that represents the wave speed in the physical domain. In other words, the pseudodifferential symbol q𝑞qitalic_q depends on the spatial variable x𝑥xitalic_x only through its dependence on the wave speed c𝑐citalic_c, which is commonly the case for models of wave propagation. The range of c𝑐citalic_c is in \mathbb{R}blackboard_R, which opens up the possibility of using univariate interpolation theory as the tool of approximation and compression. Specifically, let {c1,c2,,cM}subscript𝑐1subscript𝑐2subscript𝑐𝑀\{c_{1},c_{2},...,c_{M}\}\subset\mathbb{R}{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } ⊂ blackboard_R be a set of interpolation points and {φ1,φ2,,φM}subscript𝜑1subscript𝜑2subscript𝜑𝑀\{\varphi_{1},\varphi_{2},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_φ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } an interpolation basis of functions such that φm(cn)=0subscript𝜑𝑚subscript𝑐𝑛0\varphi_{m}(c_{n})=0italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) = 0 if nm𝑛𝑚n\neq mitalic_n ≠ italic_m and φm(cm)=1subscript𝜑𝑚subscript𝑐𝑚1\varphi_{m}(c_{m})=1italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) = 1. Examples include piecewise linear functions, Lagrange polynomials, Hermite polynomials, etc. Given this interpolation setup, we define the interpolating symbol by

qM(c,ξ)=m=1Mφm(c)q(cm,ξ),subscript𝑞𝑀𝑐𝜉superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝑐𝑞subscript𝑐𝑚𝜉\displaystyle q_{M}(c,\xi)=\sum_{m=1}^{M}\varphi_{m}(c)q(c_{m},\xi),italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ( italic_c , italic_ξ ) = ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c ) italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) , (4)

and its associated pseudodifferential operator by

(QMv)(x)=(2π)d/2m=1Mφm(c(x))q(cm,ξ)v^(ξ)eixξ𝑑ξ.subscript𝑄𝑀𝑣𝑥superscript2𝜋𝑑2superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝑐𝑥𝑞subscript𝑐𝑚𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle(Q_{M}v)(x)=(2\pi)^{-d/2}\sum_{m=1}^{M}\varphi_{m}(c(x))\int q(c_% {m},\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xi.( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c ( italic_x ) ) ∫ italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (5)

This construction allows us to approximate the action of the pseudodifferential operator (2) by the action of (5), which computationally amounts to run one fast Fourier transform (FFT) of v𝑣vitalic_v followed by M𝑀Mitalic_M inverse FFTs of the terms q(cm,ξ)v^(ξ)𝑞subscript𝑐𝑚𝜉^𝑣𝜉q(c_{m},\xi)\hat{v}(\xi)italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) for m=1,,M𝑚1𝑀m=1,...,Mitalic_m = 1 , … , italic_M. The simple, yet powerful, assumption (3) is at the heart of our proposed fast preconditioner for the Helmholtz equation. Section 2 briefly reviews some facts about pseudodifferential calculus, which we use in Section 3 to derive the pseudodifferential inversion of the Helmholtz operator and the proposed preconditioner. In Section 4 we describe the efficient evaluation of the preconditioner with nearly linear complexity with respect to the degrees of freedom and provide bounds on its accuracy based on univariate interpolation theory. In Section 5 we present some numerical experiments to illustrate the effectiveness of the preconditioner to solve the Helmholtz equation using the GMRES iterative method. We illustrate the performance of the preconditioner for a range of values for the degrees of freedom, the number M𝑀Mitalic_M of interpolation points, and the frequency of the waves. Finally, in Section 6 we provide some concluding remarks and propose some areas of potential improvement.

2 Preliminaries and notation

In this section we briefly introduce the main definitions and facts about pseudodifferential calculus. Our guiding references are [41, 27]. This section is self-contained and its notation should not be confused with notation from the rest of the paper. The first ingredient in this formulation is the Fourier transform \mathcal{F}caligraphic_F and its inverse 1superscript1\mathcal{F}^{-1}caligraphic_F start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, which for an admissible function v:d:𝑣superscript𝑑v:\mathbb{R}^{d}\to\mathbb{C}italic_v : blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT → blackboard_C, are respectively given by

v(ξ)=v^(ξ)=(2π)d/2v(x)eixξ𝑑xand1v^(x)=(2π)d/2v^(ξ)eixξ𝑑ξ.formulae-sequence𝑣𝜉^𝑣𝜉superscript2𝜋𝑑2𝑣𝑥superscript𝑒𝑖𝑥𝜉differential-d𝑥andsuperscript1^𝑣𝑥superscript2𝜋𝑑2^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle\mathcal{F}v(\xi)=\hat{v}(\xi)=(2\pi)^{-d/2}\int v(x)e^{-ix\cdot% \xi}dx\qquad\text{and}\qquad\mathcal{F}^{-1}\hat{v}(x)=(2\pi)^{-d/2}\int\hat{v% }(\xi)e^{ix\cdot\xi}d\xi.caligraphic_F italic_v ( italic_ξ ) = over^ start_ARG italic_v end_ARG ( italic_ξ ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_v ( italic_x ) italic_e start_POSTSUPERSCRIPT - italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_x and caligraphic_F start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT over^ start_ARG italic_v end_ARG ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (6)

A well-known fact about the Fourier transform is its relation with differentiation,

Dxαv(x)=(2π)d/2ξαv^(ξ)eixξ𝑑ξsubscriptsuperscript𝐷𝛼𝑥𝑣𝑥superscript2𝜋𝑑2superscript𝜉𝛼^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle D^{\alpha}_{x}v(x)=(2\pi)^{-d/2}\int\xi^{\alpha}\hat{v}(\xi)e^{% ix\cdot\xi}d\xiitalic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT italic_v ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_ξ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ (7)

for a multi-index α=(α1,,αd)d𝛼subscript𝛼1subscript𝛼𝑑superscript𝑑\alpha=(\alpha_{1},...,\alpha_{d})\in\mathbb{N}^{d}italic_α = ( italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_α start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT ) ∈ blackboard_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT where Dxα=D1α1Ddαdsubscriptsuperscript𝐷𝛼𝑥superscriptsubscript𝐷1subscript𝛼1superscriptsubscript𝐷𝑑subscript𝛼𝑑D^{\alpha}_{x}=D_{1}^{\alpha_{1}}\dots D_{d}^{\alpha_{d}}italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT = italic_D start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT … italic_D start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT end_POSTSUPERSCRIPT, Dj=ixjsubscript𝐷𝑗𝑖subscriptsubscript𝑥𝑗D_{j}=-i\partial_{x_{j}}italic_D start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT = - italic_i ∂ start_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_POSTSUBSCRIPT, and ξα=ξ1α1ξ2α2ξdαdsuperscript𝜉𝛼superscriptsubscript𝜉1subscript𝛼1superscriptsubscript𝜉2subscript𝛼2superscriptsubscript𝜉𝑑subscript𝛼𝑑\xi^{\alpha}=\xi_{1}^{\alpha_{1}}\,\xi_{2}^{\alpha_{2}}\,...\,\xi_{d}^{\alpha_% {d}}italic_ξ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT = italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT italic_ξ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT … italic_ξ start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT end_POSTSUPERSCRIPT. For a multi-index α𝛼\alphaitalic_α, we also use the notation |α|=α1++αd𝛼subscript𝛼1subscript𝛼𝑑|\alpha|=\alpha_{1}+...+\alpha_{d}| italic_α | = italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + … + italic_α start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT and α!=α1!α2!αd!𝛼subscript𝛼1subscript𝛼2subscript𝛼𝑑\alpha!=\alpha_{1}!\alpha_{2}!...\alpha_{d}!italic_α ! = italic_α start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ! italic_α start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ! … italic_α start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT !.

A differential operator of order n𝑛nitalic_n with variable coefficients qα=qα(x)subscript𝑞𝛼subscript𝑞𝛼𝑥q_{\alpha}=q_{\alpha}(x)italic_q start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT = italic_q start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ( italic_x ) can be expressed as

Q(x,D)=|α|nqα(x)Dxα.𝑄𝑥𝐷subscript𝛼𝑛subscript𝑞𝛼𝑥subscriptsuperscript𝐷𝛼𝑥\displaystyle Q(x,D)=\sum_{|\alpha|\leq n}q_{\alpha}(x)D^{\alpha}_{x}.italic_Q ( italic_x , italic_D ) = ∑ start_POSTSUBSCRIPT | italic_α | ≤ italic_n end_POSTSUBSCRIPT italic_q start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ( italic_x ) italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT . (8)

Then we have

Q(x,D)v(x)=(2π)d/2q(x,ξ)v^(ξ)eixξ𝑑ξ𝑄𝑥𝐷𝑣𝑥superscript2𝜋𝑑2𝑞𝑥𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle Q(x,D)v(x)=(2\pi)^{-d/2}\int q(x,\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xiitalic_Q ( italic_x , italic_D ) italic_v ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_q ( italic_x , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ (9)

where

q(x,ξ)=|α|nqα(x)ξα.𝑞𝑥𝜉subscript𝛼𝑛subscript𝑞𝛼𝑥superscript𝜉𝛼\displaystyle q(x,\xi)=\sum_{|\alpha|\leq n}q_{\alpha}(x)\xi^{\alpha}.italic_q ( italic_x , italic_ξ ) = ∑ start_POSTSUBSCRIPT | italic_α | ≤ italic_n end_POSTSUBSCRIPT italic_q start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ( italic_x ) italic_ξ start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT . (10)

The function q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) is called the symbol of the differential operator. The representation (9) of differential operators can be used to generalize them to a larger class known as pseudodifferential operators. Note that for a differential operator, the function q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) is polynomial with respect to ξ𝜉\xiitalic_ξ. For pseudodifferential operators, the symbols are allowed to belong to larger sets, known as Hörmander’s symbol classes defined as follows. Take n𝑛n\in\mathbb{R}italic_n ∈ blackboard_R (known as the order of the class) and define the symbol class SSnsuperscriptSS𝑛\SS^{n}roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT to consist of all Csuperscript𝐶C^{\infty}italic_C start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT functions q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) satisfying

|DxβDξαq(x,ξ)|Cαβ(1+|ξ|2)(n|α|)/2superscriptsubscript𝐷𝑥𝛽superscriptsubscript𝐷𝜉𝛼𝑞𝑥𝜉subscript𝐶𝛼𝛽superscript1superscript𝜉2𝑛𝛼2\displaystyle|D_{x}^{\beta}D_{\xi}^{\alpha}q(x,\xi)|\leq C_{\alpha\beta}(1+|% \xi|^{2})^{(n-|\alpha|)/2}| italic_D start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_β end_POSTSUPERSCRIPT italic_D start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT italic_q ( italic_x , italic_ξ ) | ≤ italic_C start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT ( italic_n - | italic_α | ) / 2 end_POSTSUPERSCRIPT (11)

for all multi-indices α𝛼\alphaitalic_α and β𝛽\betaitalic_β, and some constants Cαβ>0subscript𝐶𝛼𝛽0C_{\alpha\beta}>0italic_C start_POSTSUBSCRIPT italic_α italic_β end_POSTSUBSCRIPT > 0. For a symbol q(x,ξ)SSn𝑞𝑥𝜉superscriptSS𝑛q(x,\xi)\in\SS^{n}italic_q ( italic_x , italic_ξ ) ∈ roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT, the associated operator Op(q)Op𝑞\text{Op}\left(q\right)Op ( italic_q ), defined by

Op(q)v(x)=(2π)d/2q(x,ξ)v^(ξ)eixξ𝑑ξ,Op𝑞𝑣𝑥superscript2𝜋𝑑2𝑞𝑥𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle\text{Op}\left(q\right)v(x)=(2\pi)^{-d/2}\int q(x,\xi)\hat{v}(\xi% )e^{ix\cdot\xi}d\xi,Op ( italic_q ) italic_v ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_q ( italic_x , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ , (12)

is said to be a pseudodifferential operator that belongs to Op(SSn)OpsuperscriptSS𝑛\text{Op}\left(\SS^{n}\right)Op ( roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ). Note that by (11), when we take a derivative of a symbol q(x,ξ)𝑞𝑥𝜉q(x,\xi)italic_q ( italic_x , italic_ξ ) with respect to ξ𝜉\xiitalic_ξ, we simply obtain another pseudodifferential operator with lower order. Given a pseudodifferential operator POp(SSn)𝑃OpsuperscriptSS𝑛P\in\text{Op}\left(\SS^{n}\right)italic_P ∈ Op ( roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT ), its symbol is denoted Sym(P)Sym𝑃\text{Sym}(P)Sym ( italic_P ). Given two symbols p(x,ξ)SSn1𝑝𝑥𝜉superscriptSSsubscript𝑛1p(x,\xi)\in\SS^{n_{1}}italic_p ( italic_x , italic_ξ ) ∈ roman_SS start_POSTSUPERSCRIPT italic_n start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT and q(x,ξ)SSn2𝑞𝑥𝜉superscriptSSsubscript𝑛2q(x,\xi)\in\SS^{n_{2}}italic_q ( italic_x , italic_ξ ) ∈ roman_SS start_POSTSUPERSCRIPT italic_n start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT, the composition of their respective operators has a symbol in SSn1+n2superscriptSSsubscript𝑛1subscript𝑛2\SS^{n_{1}+n_{2}}roman_SS start_POSTSUPERSCRIPT italic_n start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_n start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUPERSCRIPT that satisfies the following relation

Sym(Op(p)Op(q))=|α|0i|α|α!Dξαp(x,ξ)Dxαq(x,ξ),mod SS,SymOp𝑝Op𝑞subscript𝛼0superscript𝑖𝛼𝛼subscriptsuperscript𝐷𝛼𝜉𝑝𝑥𝜉subscriptsuperscript𝐷𝛼𝑥𝑞𝑥𝜉mod SS\displaystyle\text{Sym}\left(\text{Op}\left(p\right)\text{Op}\left(q\right)% \right)=\sum_{|\alpha|\geq 0}\frac{i^{|\alpha|}}{\alpha!}D^{\alpha}_{\xi}p(x,% \xi)D^{\alpha}_{x}q(x,\xi),\qquad\text{mod $\SS^{-\infty}$},Sym ( Op ( italic_p ) Op ( italic_q ) ) = ∑ start_POSTSUBSCRIPT | italic_α | ≥ 0 end_POSTSUBSCRIPT divide start_ARG italic_i start_POSTSUPERSCRIPT | italic_α | end_POSTSUPERSCRIPT end_ARG start_ARG italic_α ! end_ARG italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_p ( italic_x , italic_ξ ) italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT italic_q ( italic_x , italic_ξ ) , mod roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT , (13)

where the notation mod SSsuperscriptSS\SS^{-\infty}roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT means that the difference between the left- and right-hand sides of the equality belongs to SS=nSSnsuperscriptSSsubscript𝑛superscriptSS𝑛\SS^{-\infty}=\cap_{n\in\mathbb{N}}\SS^{-n}roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT = ∩ start_POSTSUBSCRIPT italic_n ∈ blackboard_N end_POSTSUBSCRIPT roman_SS start_POSTSUPERSCRIPT - italic_n end_POSTSUPERSCRIPT. In particular, if pSSn𝑝superscriptSS𝑛p\in\SS^{n}italic_p ∈ roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT is an elliptic symbol, then p1SSnsuperscript𝑝1superscriptSS𝑛p^{-1}\in\SS^{-n}italic_p start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ∈ roman_SS start_POSTSUPERSCRIPT - italic_n end_POSTSUPERSCRIPT and

Op(p1)Op(p)=I+|α|1i|α|α!Op(Dξαp1(x,ξ)Dxαp(x,ξ)),mod SS,Opsuperscript𝑝1Op𝑝𝐼subscript𝛼1superscript𝑖𝛼𝛼Opsubscriptsuperscript𝐷𝛼𝜉superscript𝑝1𝑥𝜉subscriptsuperscript𝐷𝛼𝑥𝑝𝑥𝜉mod SS\displaystyle\text{Op}\left(p^{-1}\right)\text{Op}\left(p\right)=I+\sum_{|% \alpha|\geq 1}\frac{i^{|\alpha|}}{\alpha!}\text{Op}\left(D^{\alpha}_{\xi}p^{-1% }(x,\xi)D^{\alpha}_{x}p(x,\xi)\right),\qquad\text{mod $\SS^{-\infty}$},Op ( italic_p start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) Op ( italic_p ) = italic_I + ∑ start_POSTSUBSCRIPT | italic_α | ≥ 1 end_POSTSUBSCRIPT divide start_ARG italic_i start_POSTSUPERSCRIPT | italic_α | end_POSTSUPERSCRIPT end_ARG start_ARG italic_α ! end_ARG Op ( italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT italic_p start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ( italic_x , italic_ξ ) italic_D start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT italic_p ( italic_x , italic_ξ ) ) , mod roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT , (14)

where Dξαp1DxαpSS|α|superscriptsubscript𝐷𝜉𝛼superscript𝑝1superscriptsubscript𝐷𝑥𝛼𝑝superscriptSS𝛼D_{\xi}^{\alpha}p^{-1}D_{x}^{\alpha}p\in\SS^{-|\alpha|}italic_D start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT italic_p start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT italic_D start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT italic_p ∈ roman_SS start_POSTSUPERSCRIPT - | italic_α | end_POSTSUPERSCRIPT. Therefore, for an elliptic symbol pSSn𝑝superscriptSS𝑛p\in\SS^{n}italic_p ∈ roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT, (Op(p1)Op(p)I)SS1Opsuperscript𝑝1Op𝑝𝐼superscriptSS1(\text{Op}\left(p^{-1}\right)\text{Op}\left(p\right)-I)\in\SS^{-1}( Op ( italic_p start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) Op ( italic_p ) - italic_I ) ∈ roman_SS start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. In other words, the operator associated with the multiplicative inverse of the symbol p𝑝pitalic_p approximates the inverse of Op(p)Op𝑝\text{Op}\left(p\right)Op ( italic_p ) up to an operator of order 11-1- 1.

Also, recall that for a pseudodifferential symbol qSSn𝑞superscriptSS𝑛q\in\SS^{n}italic_q ∈ roman_SS start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT, the following Sobolev estimate

Op(q)vHs(d)2subscriptsuperscriptnormOp𝑞𝑣2superscript𝐻𝑠superscript𝑑\displaystyle\|\text{Op}\left(q\right)v\|^{2}_{H^{s}(\mathbb{R}^{d})}∥ Op ( italic_q ) italic_v ∥ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT (1+|ξ|2)ssupx|q(x,ξ)|2|v^(ξ)|2dξC002(1+|ξ|2)s+n|v^(ξ)|2𝑑ξ=C002vHs+n(d)2absentsuperscript1superscript𝜉2𝑠subscriptsupremum𝑥superscript𝑞𝑥𝜉2superscript^𝑣𝜉2𝑑𝜉superscriptsubscript𝐶002superscript1superscript𝜉2𝑠𝑛superscript^𝑣𝜉2differential-d𝜉superscriptsubscript𝐶002subscriptsuperscriptnorm𝑣2superscript𝐻𝑠𝑛superscript𝑑\displaystyle\leq\int(1+|\xi|^{2})^{s}\sup_{x}{|q(x,\xi)|^{2}}|\hat{v}(\xi)|^{% 2}d\xi\leq C_{00}^{2}\int(1+|\xi|^{2})^{s+n}|\hat{v}(\xi)|^{2}d\xi=C_{00}^{2}% \|v\|^{2}_{H^{s+n}(\mathbb{R}^{d})}≤ ∫ ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT roman_sup start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT | italic_q ( italic_x , italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | over^ start_ARG italic_v end_ARG ( italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d italic_ξ ≤ italic_C start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∫ ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT | over^ start_ARG italic_v end_ARG ( italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d italic_ξ = italic_C start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∥ italic_v ∥ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT (15)

holds provided that vHs+n(d)𝑣superscript𝐻𝑠𝑛superscript𝑑v\in H^{s+n}(\mathbb{R}^{d})italic_v ∈ italic_H start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) and where C00>0subscript𝐶000C_{00}>0italic_C start_POSTSUBSCRIPT 00 end_POSTSUBSCRIPT > 0 is the constant appearing in (11). In other words, Op(q)Op𝑞\text{Op}\left(q\right)Op ( italic_q ) maps Hs+n(d)superscript𝐻𝑠𝑛superscript𝑑H^{s+n}(\mathbb{R}^{d})italic_H start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) continuously into Hs(d)superscript𝐻𝑠superscript𝑑H^{s}(\mathbb{R}^{d})italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) for any s𝑠s\in\mathbb{R}italic_s ∈ blackboard_R. Moreover, the norm of Op(q)Op𝑞\text{Op}\left(q\right)Op ( italic_q ) is bounded by the constant appearing in (11) that bounds its symbol q𝑞qitalic_q.

3 Pseudodifferential inversion of the Helmholtz equation

We are interested in efficiently computing approximations to the solution u𝑢uitalic_u of the Helmholtz equation (1) subject to the Sommerfeld radiation condition. In practice, we can only approximately compute the solution u𝑢uitalic_u restricted to a bounded domain ΩdΩsuperscript𝑑\Omega\subset\mathbb{R}^{d}roman_Ω ⊂ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT with Lipschitz boundary ΩΩ\partial\Omega∂ roman_Ω. In the Helmholtz equation (1), Δ=x12++xd2Δsuperscriptsubscriptsubscript𝑥12superscriptsubscriptsubscript𝑥𝑑2\Delta=\partial_{x_{1}}^{2}+...+\partial_{x_{d}}^{2}roman_Δ = ∂ start_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + … + ∂ start_POSTSUBSCRIPT italic_x start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT denotes the Laplacian, ω𝜔\omegaitalic_ω is the angular frequency of oscillation, c=c(x)𝑐𝑐𝑥c=c(x)italic_c = italic_c ( italic_x ) is the variable wave speed, a=a(x)𝑎𝑎𝑥a=a(x)italic_a = italic_a ( italic_x ) is the damping coefficient, and f𝑓fitalic_f is a given source with support in ΩΩ\Omegaroman_Ω. We assume that the wave speed c(x)𝑐𝑥c(x)italic_c ( italic_x ) is smooth, and bounded from above and below, that is, 0<cminccmax<0subscript𝑐min𝑐subscript𝑐max0<c_{\rm min}\leq c\leq c_{\rm max}<\infty0 < italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT ≤ italic_c ≤ italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < ∞. Similarly, we assume that the damping coefficient a(x)𝑎𝑥a(x)italic_a ( italic_x ) is a smooth function such that 0<aminaamax<0subscript𝑎min𝑎subscript𝑎max0<a_{\rm min}\leq a\leq a_{\rm max}<\infty0 < italic_a start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT ≤ italic_a ≤ italic_a start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < ∞. We also assume that there are constant background properties cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT and aosubscript𝑎𝑜a_{o}italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT such that (c(x)co)𝑐𝑥subscript𝑐𝑜(c(x)-c_{o})( italic_c ( italic_x ) - italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ) and (a(x)ao)𝑎𝑥subscript𝑎𝑜(a(x)-a_{o})( italic_a ( italic_x ) - italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ) are compactly supported in ΩΩ\Omegaroman_Ω. The well-posedness of this problem in both weak and strong formulations has been established. See for instance [33, 11].

Our approach relies on using pseudodifferential calculus to approximate the inverse of the Helmholtz operator governing (1). The Helmholtz operator is defined by

P=ω2+iωa+c2ΔOp(SS2),𝑃superscript𝜔2𝑖𝜔𝑎superscript𝑐2ΔOpsuperscriptSS2\displaystyle P=\omega^{2}+i\omega a+c^{2}\Delta\in\text{Op}\left(\SS^{2}% \right),italic_P = italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a + italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ ∈ Op ( roman_SS start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) , (16)

and its pseudodifferential symbol is given by

Sym(P)=p=ω2+iωac2(ξ12++ξd2)=ω2+iωac2|ξ|2SS2,Sym𝑃𝑝superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscriptsubscript𝜉12superscriptsubscript𝜉𝑑2superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉2superscriptSS2\displaystyle\text{Sym}\left(P\right)=p=\omega^{2}+i\omega a-c^{2}\left(\xi_{1% }^{2}+...+\xi_{d}^{2}\right)=\omega^{2}+i\omega a-c^{2}|\xi|^{2}\in\SS^{2},Sym ( italic_P ) = italic_p = italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + … + italic_ξ start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) = italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∈ roman_SS start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT , (17)

where the Fourier dual of the spatial variable x=(x1,,xd)d𝑥subscript𝑥1subscript𝑥𝑑superscript𝑑x=(x_{1},...,x_{d})\in\mathbb{R}^{d}italic_x = ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_x start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT ) ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT is denoted by ξ=(ξ1,,ξd)d𝜉subscript𝜉1subscript𝜉𝑑superscript𝑑\xi=(\xi_{1},...,\xi_{d})\in\mathbb{R}^{d}italic_ξ = ( italic_ξ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_ξ start_POSTSUBSCRIPT italic_d end_POSTSUBSCRIPT ) ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT. Due to the above assumptions on the wave speed and damping, it can be shown that the operator P1Op(SS2)superscript𝑃1OpsuperscriptSS2P^{-1}\in\text{Op}\left(\SS^{-2}\right)italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ∈ Op ( roman_SS start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT ) exists [41, 27]. Our objective is to construct an approximation for the symbol of P1superscript𝑃1P^{-1}italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. First, using the formula (13) for the symbol of the composition of the operators POp(SS2)𝑃OpsuperscriptSS2P\in\text{Op}\left(\SS^{2}\right)italic_P ∈ Op ( roman_SS start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) and P1Op(SS2)superscript𝑃1OpsuperscriptSS2P^{-1}\in\text{Op}\left(\SS^{-2}\right)italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ∈ Op ( roman_SS start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT ), we obtain

1=Sym(PP1)=|α|0i|α|α!DξαSym(P)DxαSym(P1),mod SS.formulae-sequence1Sym𝑃superscript𝑃1subscript𝛼0superscript𝑖𝛼𝛼superscriptsubscript𝐷𝜉𝛼Sym𝑃superscriptsubscript𝐷𝑥𝛼Symsuperscript𝑃1mod SS\displaystyle 1=\text{Sym}\left(PP^{-1}\right)=\sum_{|\alpha|\geq 0}\frac{i^{|% \alpha|}}{\alpha!}D_{\xi}^{\alpha}\text{Sym}\left(P\right)D_{x}^{\alpha}\text{% Sym}\left(P^{-1}\right),\qquad\text{mod $\SS^{-\infty}$}.1 = Sym ( italic_P italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) = ∑ start_POSTSUBSCRIPT | italic_α | ≥ 0 end_POSTSUBSCRIPT divide start_ARG italic_i start_POSTSUPERSCRIPT | italic_α | end_POSTSUPERSCRIPT end_ARG start_ARG italic_α ! end_ARG italic_D start_POSTSUBSCRIPT italic_ξ end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT Sym ( italic_P ) italic_D start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_α end_POSTSUPERSCRIPT Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) , mod roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT . (18)

Note that the expansion on the right-hand side of (18) contains only 3 terms because Sym(P)Sym𝑃\text{Sym}\left(P\right)Sym ( italic_P ), given by (17), is a second degree polynomial with respect to ξ𝜉\xiitalic_ξ. This implies that Sym(P1)Symsuperscript𝑃1\text{Sym}\left(P^{-1}\right)Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) satisfies the following second-order linear partial differential equation

c2ΔSym(P1)+2ic2ξSym(P1)+(ω2+iωac2|ξ|2)Sym(P1)=1,mod SS,superscript𝑐2ΔSymsuperscript𝑃12𝑖superscript𝑐2𝜉Symsuperscript𝑃1superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉2Symsuperscript𝑃11mod SS\displaystyle c^{2}\Delta\,\text{Sym}\left(P^{-1}\right)+2ic^{2}\xi\cdot\nabla% \text{Sym}\left(P^{-1}\right)+\left(\omega^{2}+i\omega a-c^{2}|\xi|^{2}\right)% \,\text{Sym}\left(P^{-1}\right)=1,\qquad\text{mod $\SS^{-\infty}$},italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) + 2 italic_i italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ ⋅ ∇ Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) + ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) = 1 , mod roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT , (19)

as a function of the space variables xd𝑥superscript𝑑x\in\mathbb{R}^{d}italic_x ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT for each value of ξd𝜉superscript𝑑\xi\in\mathbb{R}^{d}italic_ξ ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT.

Now we propose the following expansion for the symbol of P1superscript𝑃1P^{-1}italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT,

q=j=2+qj,mod SS𝑞superscriptsubscript𝑗2subscript𝑞𝑗mod SS\displaystyle q=\sum_{j=2}^{+\infty}q_{-j},\qquad\text{mod $\SS^{-\infty}$}italic_q = ∑ start_POSTSUBSCRIPT italic_j = 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT + ∞ end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT - italic_j end_POSTSUBSCRIPT , mod roman_SS start_POSTSUPERSCRIPT - ∞ end_POSTSUPERSCRIPT (20)

where qjSSjsubscript𝑞𝑗superscriptSS𝑗q_{-j}\in\SS^{-j}italic_q start_POSTSUBSCRIPT - italic_j end_POSTSUBSCRIPT ∈ roman_SS start_POSTSUPERSCRIPT - italic_j end_POSTSUPERSCRIPT and where the principal symbol is defined as

q2=1p=1ω2+iωac2|ξ|2.subscript𝑞21𝑝1superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉2\displaystyle q_{-2}=\frac{1}{p}=\frac{1}{\omega^{2}+i\omega a-c^{2}|\xi|^{2}}.italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT = divide start_ARG 1 end_ARG start_ARG italic_p end_ARG = divide start_ARG 1 end_ARG start_ARG italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (21)

Plugging (20)-(21) into (19) and gathering terms with the same classes, we define all of the symbols recursively

q3subscript𝑞3\displaystyle q_{-3}italic_q start_POSTSUBSCRIPT - 3 end_POSTSUBSCRIPT =2ic2ξq2ω2+iωac2|ξ|2=2ic2ξ(iωa|ξ|2c2)(ω2+iωac2|ξ|2)3,absent2𝑖superscript𝑐2𝜉subscript𝑞2superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉22𝑖superscript𝑐2𝜉𝑖𝜔𝑎superscript𝜉2superscript𝑐2superscriptsuperscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉23\displaystyle=-\frac{2ic^{2}\xi\cdot\nabla q_{-2}}{\omega^{2}+i\omega a-c^{2}|% \xi|^{2}}=\frac{2ic^{2}\xi\cdot\left(i\omega\nabla a-|\xi|^{2}\nabla c^{2}% \right)}{\left(\omega^{2}+i\omega a-c^{2}|\xi|^{2}\right)^{3}},= - divide start_ARG 2 italic_i italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ ⋅ ∇ italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG = divide start_ARG 2 italic_i italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ ⋅ ( italic_i italic_ω ∇ italic_a - | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∇ italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG start_ARG ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG , (22)
qjsubscript𝑞𝑗\displaystyle q_{-j}italic_q start_POSTSUBSCRIPT - italic_j end_POSTSUBSCRIPT =2ic2ξqj+1+c2Δqj+2ω2+iωac2|ξ|2,j4.formulae-sequenceabsent2𝑖superscript𝑐2𝜉subscript𝑞𝑗1superscript𝑐2Δsubscript𝑞𝑗2superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉2𝑗4\displaystyle=-\frac{2ic^{2}\xi\cdot\nabla q_{-j+1}+c^{2}\Delta q_{-j+2}}{% \omega^{2}+i\omega a-c^{2}|\xi|^{2}},\qquad j\geq 4.= - divide start_ARG 2 italic_i italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ξ ⋅ ∇ italic_q start_POSTSUBSCRIPT - italic_j + 1 end_POSTSUBSCRIPT + italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT roman_Δ italic_q start_POSTSUBSCRIPT - italic_j + 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_j ≥ 4 . (23)

As opposed to a classical expansion of the symbol Sym(P1)Symsuperscript𝑃1\text{Sym}\left(P^{-1}\right)Sym ( italic_P start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ), the adhoc expansion (20) contains exactly one non-vanishing term when the acoustic medium is homogeneous. Hence, we expect a truncation of this expansion to be more accurate than a truncated classical expansion when the medium is homogeneous or when there are small variations in the wave speed c𝑐citalic_c and damping a𝑎aitalic_a.

The proposed preconditioner is the pseudodifferential operator associated with the principal symbol of the expansion (20),

(Qv)(x)=(2π)d/2q2(x,ξ)v^(ξ)eixξ𝑑ξ.𝑄𝑣𝑥superscript2𝜋𝑑2subscript𝑞2𝑥𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle(Qv)(x)=(2\pi)^{-d/2}\int q_{-2}(x,\xi)\hat{v}(\xi)e^{ix\cdot\xi}% d\xi.( italic_Q italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT ( italic_x , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (24)

Note that the symbol pSS2𝑝superscriptSS2p\in\SS^{2}italic_p ∈ roman_SS start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT defined in (17) is clearly elliptic [41, 27]. Therefore, following (14), we obtain that QPI=ROp(SS1)𝑄𝑃𝐼𝑅OpsuperscriptSS1QP-I=R\in\text{Op}\left(\SS^{-1}\right)italic_Q italic_P - italic_I = italic_R ∈ Op ( roman_SS start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ). This means that the proposed preconditioner Q𝑄Qitalic_Q defined by (24) approximates the inverse of the Helmholtz operator up to an error R𝑅Ritalic_R with pseudodifferential order of 11-1- 1. By definition, such an operator possesses a symbol r𝑟ritalic_r that decays in the frequency domain, ie., |r|C(1+|ξ|2)1/2C/|ξ|𝑟𝐶superscript1superscript𝜉212similar-to𝐶𝜉|r|\leq C(1+|\xi|^{2})^{-1/2}\sim C/|\xi|| italic_r | ≤ italic_C ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT ∼ italic_C / | italic_ξ | as |ξ|𝜉|\xi|\to\infty| italic_ξ | → ∞. Consequently, we would expect the preconditioner to be more effective when restricted to problems whose solutions are highly oscillatory. We observe numerical evidence of this behavior in Section 5 below.

We may also characterize the performance of the preconditioner in terms of its spectral behavior. For this purpose, we consider the solutions u𝑢uitalic_u of (1) restricted to the bounded domain ΩΩ\Omegaroman_Ω using the restriction operator ΩsubscriptΩ\mathcal{E}_{\Omega}caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT that maps Hs(d)uu|ΩHs(Ω)containssuperscript𝐻𝑠superscript𝑑𝑢maps-toevaluated-at𝑢Ωsuperscript𝐻𝑠ΩH^{s}(\mathbb{R}^{d})\owns u\mapsto u|_{\Omega}\in H^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) ∋ italic_u ↦ italic_u | start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ∈ italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ). Now note that since u𝑢uitalic_u is a radiating solution to the homogeneous Helmholtz equation in dΩsuperscript𝑑Ω\mathbb{R}^{d}\setminus\Omegablackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ∖ roman_Ω (because f𝑓fitalic_f is supported in ΩΩ\Omegaroman_Ω), then Pu𝑃𝑢Puitalic_P italic_u has support in ΩΩ\Omegaroman_Ω and u|dΩevaluated-at𝑢superscript𝑑Ωu|_{\mathbb{R}^{d}\setminus\Omega}italic_u | start_POSTSUBSCRIPT blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ∖ roman_Ω end_POSTSUBSCRIPT is fully determined by u|Ωevaluated-at𝑢Ωu|_{\Omega}italic_u | start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT. Therefore, we can think of P𝑃Pitalic_P as mapping Hs(Ω)superscript𝐻𝑠ΩH^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) into Hs2(Ω)superscript𝐻𝑠2ΩH^{s-2}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s - 2 end_POSTSUPERSCRIPT ( roman_Ω ). The operator RΩ=ΩQPI:Hs(Ω)Hs(Ω):subscript𝑅ΩsubscriptΩ𝑄𝑃𝐼superscript𝐻𝑠Ωsuperscript𝐻𝑠ΩR_{\Omega}=\mathcal{E}_{\Omega}QP-I:H^{s}(\Omega)\to H^{s}(\Omega)italic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT = caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P - italic_I : italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) still has pseudodifferential order of 11-1- 1. See details in [27, Ch. 4]. Hence, for any Sobolev scale s𝑠s\in\mathbb{R}italic_s ∈ blackboard_R, the operator RΩ=ΩQPIsubscript𝑅ΩsubscriptΩ𝑄𝑃𝐼R_{\Omega}=\mathcal{E}_{\Omega}QP-Iitalic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT = caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P - italic_I is bounded from Hs(Ω)superscript𝐻𝑠ΩH^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) to Hs+1(Ω)superscript𝐻𝑠1ΩH^{s+1}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s + 1 end_POSTSUPERSCRIPT ( roman_Ω ) and compact from Hs(Ω)superscript𝐻𝑠ΩH^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) to Hs(Ω)superscript𝐻𝑠ΩH^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) due to the compact embedding Hs+1(Ω)Hs(Ω)double-subset-ofsuperscript𝐻𝑠1Ωsuperscript𝐻𝑠ΩH^{s+1}(\Omega)\Subset H^{s}(\Omega)italic_H start_POSTSUPERSCRIPT italic_s + 1 end_POSTSUPERSCRIPT ( roman_Ω ) ⋐ italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) for a bounded domain ΩΩ\Omegaroman_Ω [33, Thm. 3.27]. As a direct consequence of the spectral behavior of compact operators ([32, Thm 3.9], [14, Corollary 2.2.13], [30, Thm 6.26, Ch 3]), we can make the following assertion.

Theorem 1.

The preconditioned operator ΩQP:Hs(Ω)Hs(Ω):subscriptΩ𝑄𝑃superscript𝐻𝑠Ωsuperscript𝐻𝑠Ω\mathcal{E}_{\Omega}QP:H^{s}(\Omega)\to H^{s}(\Omega)caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P : italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) is Fredholm (being a compact perturbation of the identity), and its spectrum consists of z=1𝑧1z=1italic_z = 1 and at most a countable set of eigenvalues with no point of accumulation except, possibly, z=1𝑧1z=1italic_z = 1.

This theorem characterizes the spectral performance of the preconditioner Q𝑄Qitalic_Q at the continuous level. In computational practice, both Q𝑄Qitalic_Q and P𝑃Pitalic_P are discretized and it is not clear at this point how Theorem 1 translates into the discrete setting. However, Theorem 1 provides an insight regarding the ability of the preconditioner Q𝑄Qitalic_Q to cluster the eigenvalues of QP𝑄𝑃QPitalic_Q italic_P, away from zero and infinity, in order to improve the convergence of the GMRES iterative method [42, Part VI].

4 Evaluation algorithm

Now we describe the proposed algorithm to evaluate the pseudodifferential preconditioner defined by (24), its expected computational complexity and accuracy. We propose an algorithm for a pseudodifferential operator of the following form

(Qv)(x)=(2π)d/2g(x)q(c(x),ξ)v^(ξ)eixξ𝑑ξ,𝑄𝑣𝑥superscript2𝜋𝑑2𝑔𝑥𝑞𝑐𝑥𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle(Qv)(x)=(2\pi)^{-d/2}g(x)\int q(c(x),\xi)\hat{v}(\xi)e^{ix\cdot% \xi}d\xi,( italic_Q italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT italic_g ( italic_x ) ∫ italic_q ( italic_c ( italic_x ) , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ , (25)

where the symbol q𝑞qitalic_q satisfies (3), that is, it depends on the spatial variable xd𝑥superscript𝑑x\in\mathbb{R}^{d}italic_x ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT through the function c=c(x)𝑐𝑐𝑥c=c(x)italic_c = italic_c ( italic_x ) which plays the role of the heterogeneous wave speed in acoustic problems. For the principal symbol (21), this assumption requires that the damping coefficient a=a(c(x))𝑎𝑎𝑐𝑥a=a(c(x))italic_a = italic_a ( italic_c ( italic_x ) ) be a function of the wave speed c𝑐citalic_c. This means that physical media with the same wave speed have the same damping coefficient as well. Also note that multiplication by g(x)𝑔𝑥g(x)italic_g ( italic_x ) in (25) does not pose any additional computational difficulties. Hence, in the sequel we set g(x)=1𝑔𝑥1g(x)=1italic_g ( italic_x ) = 1 for simplicity.

4.1 Algorithm

As stated in the previous section, it is assumed that c=c(x)𝑐𝑐𝑥c=c(x)italic_c = italic_c ( italic_x ) is a real-valued smooth function bounded from below and from above such that 0<cminccmax<0subscript𝑐min𝑐subscript𝑐max0<c_{\rm min}\leq c\leq c_{\rm max}<\infty0 < italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT ≤ italic_c ≤ italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT < ∞. Let {c1,c2,,cM}[cmin,cmax]subscript𝑐1subscript𝑐2subscript𝑐𝑀subscript𝑐minsubscript𝑐max\{c_{1},c_{2},...,c_{M}\}\subset[c_{\rm min},c_{\rm max}]{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } ⊂ [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] be a finite ordered set of distinct interpolation points such that c1=cminsubscript𝑐1subscript𝑐minc_{1}=c_{\rm min}italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT and cM=cmaxsubscript𝑐𝑀subscript𝑐maxc_{M}=c_{\rm max}italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT. Also let {φ1,,φM}subscript𝜑1subscript𝜑𝑀\{\varphi_{1},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } be an interpolating basis associated with the points {c1,c2,,cM}subscript𝑐1subscript𝑐2subscript𝑐𝑀\{c_{1},c_{2},...,c_{M}\}{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT }, that is, each φm:[cmin,cmax][0,1]:subscript𝜑𝑚subscript𝑐minsubscript𝑐max01\varphi_{m}:[c_{\rm min},c_{\rm max}]\to[0,1]italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT : [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] → [ 0 , 1 ] is a function such that φm(cn)=0subscript𝜑𝑚subscript𝑐𝑛0\varphi_{m}(c_{n})=0italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) = 0 for all nm𝑛𝑚n\neq mitalic_n ≠ italic_m, and φm(cm)=1subscript𝜑𝑚subscript𝑐𝑚1\varphi_{m}(c_{m})=1italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) = 1 for each m=1,,M𝑚1𝑀m=1,...,Mitalic_m = 1 , … , italic_M. Using this interpolating basis, we define the pseudodifferential interpolating symbol as

qM(c,ξ)=m=1Mφm(c)q(cm,ξ),subscript𝑞𝑀𝑐𝜉superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝑐𝑞subscript𝑐𝑚𝜉\displaystyle q_{M}(c,\xi)=\sum_{m=1}^{M}\varphi_{m}(c)q(c_{m},\xi),italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ( italic_c , italic_ξ ) = ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c ) italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) , (26)

and its associated pseudodifferential operator as

(QMv)(x)=(2π)d/2m=1Mφm(c(x))q(cm,ξ)v^(ξ)eixξ𝑑ξ.subscript𝑄𝑀𝑣𝑥superscript2𝜋𝑑2superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝑐𝑥𝑞subscript𝑐𝑚𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle(Q_{M}v)(x)=(2\pi)^{-d/2}\sum_{m=1}^{M}\varphi_{m}(c(x))\int q(c_% {m},\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xi.( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c ( italic_x ) ) ∫ italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (27)

4.2 Complexity estimates

The complexity of the algorithm will be analyzed in terms of the number of multiplications and the discretization of the spatial domain. For simplicity, let us consider a regular d𝑑ditalic_d-dimensional grid of a hybercube containing the support of the input function v=v(x)𝑣𝑣𝑥v=v(x)italic_v = italic_v ( italic_x ). Let N𝑁Nitalic_N denote the number of grid points in each coordinate direction. Hence, there are Ndsuperscript𝑁𝑑N^{d}italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT discrete evaluation points known as degrees of freedom (DOF). We assume that the discrete Fourier transform and the inverse discrete Fourier transform are computed using FFT-type algorithms with one-dimensional complexity 𝒪(NlogN)𝒪𝑁𝑁\mathcal{O}\left(N\log N\right)caligraphic_O ( italic_N roman_log italic_N ). Hence, we expect the direct computational evaluation of a pseudodifferential operator of the form (25) to have 𝒪(N2dlogN)𝒪superscript𝑁2𝑑𝑁\mathcal{O}\left(N^{2d}\log N\right)caligraphic_O ( italic_N start_POSTSUPERSCRIPT 2 italic_d end_POSTSUPERSCRIPT roman_log italic_N ) complexity.

For the proposed algorithm, the practical evaluation of (27) amounts to compute one FFT of v(x)𝑣𝑥v(x)italic_v ( italic_x ) followed by the sum of M𝑀Mitalic_M inverse FFTs, one for each input q(cm,ξ)v^(ξ)𝑞subscript𝑐𝑚𝜉^𝑣𝜉q(c_{m},\xi)\hat{v}(\xi)italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) for m=1,,M𝑚1𝑀m=1,...,Mitalic_m = 1 , … , italic_M. Hence, we obtain the following result.

Lemma 1.

The proposed algorithm based on (27) to approximate the action of a pseudodifferential operator of the form (25) has 𝒪(MNdlogN)𝒪𝑀superscript𝑁𝑑𝑁\mathcal{O}\left(MN^{d}\log N\right)caligraphic_O ( italic_M italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT roman_log italic_N ) complexity.

Since DOF=Ndsuperscript𝑁𝑑N^{d}italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT, the proposed algorithm is nearly linear with respect to the degrees of freedom. Hence, in general, we expect the proposed algorithm to be much more computationally efficient than the naive algorithm when MNdmuch-less-than𝑀superscript𝑁𝑑M\ll N^{d}italic_M ≪ italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT, provided that M𝑀Mitalic_M is sufficiently large to obtain satisfactory accuracy for the interpolation of the pseudodifferential symbol q𝑞qitalic_q. The theory of interpolation dictates that, for a well-chosen interpolation basis {φ1,,φM}subscript𝜑1subscript𝜑𝑀\{\varphi_{1},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT }, the number M𝑀Mitalic_M depends only on the smoothness of the symbol q𝑞qitalic_q to render accuracy in the sup norm. Hence, when q𝑞qitalic_q is sufficiently smooth, we can expect M𝑀Mitalic_M to be independent of N𝑁Nitalic_N. This is extremely advantageous in situations where v𝑣vitalic_v is highly oscillatory or when q𝑞qitalic_q promotes high frequencies, thus requiring N𝑁Nitalic_N to be large to resolve such oscillations. For instance, in a typical setting for ultrasound propagation in biomedical applications, there are about 100similar-toabsent100\sim 100∼ 100 wavelengths across the domain. If each wavelength is discretized with about 10similar-toabsent10\sim 10∼ 10 points, then N1000similar-to𝑁1000N\sim 1000italic_N ∼ 1000. Therefore, DOF=Nd106similar-tosuperscript𝑁𝑑superscript106N^{d}\sim 10^{6}italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ∼ 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT in a two-dimensional setting (d=2𝑑2d=2italic_d = 2). This number of DOF is many orders of magnitude greater than the number of interpolation points M10similar-to𝑀10M\sim 10italic_M ∼ 10 needed to resolve the profile of a pseudodifferential symbol like (21) with reasonable accuracy. We characterize the accuracy of (27) in more detail as follows.

4.3 Accuracy and conditioning estimates

The accuracy of the proposed algorithm (26)-(27) depends on its ability to approximate the pseudodifferential symbol q𝑞qitalic_q of the operator (25). A bound on the error |qqM|𝑞subscript𝑞𝑀|q-q_{M}|| italic_q - italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT | can be established depending on the choice of interpolation points {c1,c2,,cM}subscript𝑐1subscript𝑐2subscript𝑐𝑀\{c_{1},c_{2},...,c_{M}\}{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } and interpolation functions {φ1,φ2,,φM}subscript𝜑1subscript𝜑2subscript𝜑𝑀\{\varphi_{1},\varphi_{2},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_φ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT }. Here we describe error bounds for the piecewise linear interpolation strategy. For details and other polynomial strategies, see [31, Ch. 8] or [5, Ch. 3]. The proof follows standard arguments for the real and imaginary parts of the difference (qqM)𝑞subscript𝑞𝑀(q-q_{M})( italic_q - italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ).

Lemma 2.

Let q:[cmin,cmax]×d:𝑞subscript𝑐minsubscript𝑐maxsuperscript𝑑q:[c_{\rm min},c_{\rm max}]\times\mathbb{R}^{d}\to\mathbb{C}italic_q : [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] × blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT → blackboard_C be a smooth function such that

|2qc2(c,ξ)|C(1+|ξ|2)n/2,superscript2𝑞superscript𝑐2𝑐𝜉𝐶superscript1superscript𝜉2𝑛2\displaystyle\left|\frac{\partial^{2}q}{\partial c^{2}}(c,\xi)\right|\leq C(1+% |\xi|^{2})^{n/2},| divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q end_ARG start_ARG ∂ italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ( italic_c , italic_ξ ) | ≤ italic_C ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_n / 2 end_POSTSUPERSCRIPT , (28)

for some n𝑛n\in\mathbb{R}italic_n ∈ blackboard_R, and for some constant C>0𝐶0C>0italic_C > 0 independent of ξd𝜉superscript𝑑\xi\in\mathbb{R}^{d}italic_ξ ∈ blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT. If {φ1,φ2,,φM}subscript𝜑1subscript𝜑2subscript𝜑𝑀\{\varphi_{1},\varphi_{2},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_φ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } is a piecewise linear interpolation basis for the equidistant interpolation points {c1,c2,,cM}subscript𝑐1subscript𝑐2subscript𝑐𝑀\{c_{1},c_{2},...,c_{M}\}{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } with step size hM=(cmaxcmin)/(M1)subscript𝑀subscript𝑐maxsubscript𝑐min𝑀1h_{M}=(c_{\rm max}-c_{\rm min})/(M-1)italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT = ( italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT - italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT ) / ( italic_M - 1 ), then

|q(c,ξ)qM(c,ξ)|ChM242(1+|ξ|2)n/2.𝑞𝑐𝜉subscript𝑞𝑀𝑐𝜉𝐶superscriptsubscript𝑀242superscript1superscript𝜉2𝑛2\displaystyle\left|q(c,\xi)-q_{M}(c,\xi)\right|\leq\frac{Ch_{M}^{2}}{4\sqrt{2}% }(1+|\xi|^{2})^{n/2}.| italic_q ( italic_c , italic_ξ ) - italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ( italic_c , italic_ξ ) | ≤ divide start_ARG italic_C italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 square-root start_ARG 2 end_ARG end_ARG ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_n / 2 end_POSTSUPERSCRIPT . (29)

As a result of this estimate, we can easily obtain the following error bound for the evaluation of the pseudodifferential operators (25)-(27).

Theorem 2.

Let the symbol q𝑞qitalic_q, {φ1,φ2,,φM}subscript𝜑1subscript𝜑2subscript𝜑𝑀\{\varphi_{1},\varphi_{2},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_φ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } and hMsubscript𝑀h_{M}italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT satisfy the assumptions of Lemma 2, and let s𝑠s\in\mathbb{R}italic_s ∈ blackboard_R denote a Sobolev scale. Then

(QQM)vHs(d)ChM242vHs+n(d).subscriptnorm𝑄subscript𝑄𝑀𝑣superscript𝐻𝑠superscript𝑑𝐶superscriptsubscript𝑀242subscriptnorm𝑣superscript𝐻𝑠𝑛superscript𝑑\displaystyle\|(Q-Q_{M})v\|_{H^{s}(\mathbb{R}^{d})}\leq\frac{Ch_{M}^{2}}{4% \sqrt{2}}\|v\|_{H^{s+n}(\mathbb{R}^{d})}.∥ ( italic_Q - italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ) italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT ≤ divide start_ARG italic_C italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 square-root start_ARG 2 end_ARG end_ARG ∥ italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT . (30)
Proof.

To obtain the estimate (30), we simply employ the definition of the pseudodifferential operators, the equivalence of Sobolev norms in the Fourier space, and the estimate (29) from Lemma 2 to obtain,

(QQM)vHs(d)2superscriptsubscriptnorm𝑄subscript𝑄𝑀𝑣superscript𝐻𝑠superscript𝑑2\displaystyle\|(Q-Q_{M})v\|_{H^{s}(\mathbb{R}^{d})}^{2}∥ ( italic_Q - italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ) italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (1+|ξ|2)ssupx|q(c(x),ξ)qM(c(x),ξ)|2|v^(ξ)|2dξabsentsuperscript1superscript𝜉2𝑠subscriptsupremum𝑥superscript𝑞𝑐𝑥𝜉subscript𝑞𝑀𝑐𝑥𝜉2superscript^𝑣𝜉2𝑑𝜉\displaystyle\leq\int(1+|\xi|^{2})^{s}\sup_{x}\left|q(c(x),\xi)-q_{M}(c(x),\xi% )\right|^{2}|\hat{v}(\xi)|^{2}d\xi≤ ∫ ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT roman_sup start_POSTSUBSCRIPT italic_x end_POSTSUBSCRIPT | italic_q ( italic_c ( italic_x ) , italic_ξ ) - italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ( italic_c ( italic_x ) , italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | over^ start_ARG italic_v end_ARG ( italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d italic_ξ
(ChM242)2(1+|ξ|2)s+n|v^(ξ)|2𝑑ξabsentsuperscript𝐶superscriptsubscript𝑀2422superscript1superscript𝜉2𝑠𝑛superscript^𝑣𝜉2differential-d𝜉\displaystyle\leq\left(\frac{Ch_{M}^{2}}{4\sqrt{2}}\right)^{2}\int(1+|\xi|^{2}% )^{s+n}|\hat{v}(\xi)|^{2}d\xi≤ ( divide start_ARG italic_C italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 square-root start_ARG 2 end_ARG end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∫ ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT | over^ start_ARG italic_v end_ARG ( italic_ξ ) | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_d italic_ξ
(ChM242)2vHs+n(d)2.absentsuperscript𝐶superscriptsubscript𝑀2422superscriptsubscriptnorm𝑣superscript𝐻𝑠𝑛superscript𝑑2\displaystyle\leq\left(\frac{Ch_{M}^{2}}{4\sqrt{2}}\right)^{2}\|v\|_{H^{s+n}(% \mathbb{R}^{d})}^{2}.≤ ( divide start_ARG italic_C italic_h start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 4 square-root start_ARG 2 end_ARG end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ∥ italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s + italic_n end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT .

which renders the desired result for the same constant C>0𝐶0C>0italic_C > 0 appearing in the assumption (28). ∎

Finally, we conclude this section by characterizing the ability of the proposed preconditioner QMsubscript𝑄𝑀Q_{M}italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT to approximately invert the Helmholtz operator P𝑃Pitalic_P. First, we verify that the assumptions of Lemma 2 are indeed satisfied by the principal symbol (21). By assumption, the wave speed c𝑐citalic_c is smooth, bounded above and below, and (cco)𝑐subscript𝑐o(c-c_{\rm o})( italic_c - italic_c start_POSTSUBSCRIPT roman_o end_POSTSUBSCRIPT ) has compact support. So c𝑐citalic_c attains a minimum and maximum which we denote by cminsubscript𝑐minc_{\rm min}italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT and cmaxsubscript𝑐maxc_{\rm max}italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, respectively. The derivatives of (21) with respect to c𝑐citalic_c are as follows

q2csubscript𝑞2𝑐\displaystyle\frac{\partial q_{-2}}{\partial c}divide start_ARG ∂ italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_c end_ARG =2c|ξ|2(ω2+iωac2|ξ|2)2,absent2𝑐superscript𝜉2superscriptsuperscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉22\displaystyle=\frac{2c|\xi|^{2}}{\left(\omega^{2}+i\omega a-c^{2}|\xi|^{2}% \right)^{2}},= divide start_ARG 2 italic_c | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ,
2q2c2superscript2subscript𝑞2superscript𝑐2\displaystyle\frac{\partial^{2}q_{-2}}{\partial c^{2}}divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT end_ARG start_ARG ∂ italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG =2|ξ|2(ω2+iωac2|ξ|2)+8c2|ξ|4(ω2+iωac2|ξ|2)3.absent2superscript𝜉2superscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉28superscript𝑐2superscript𝜉4superscriptsuperscript𝜔2𝑖𝜔𝑎superscript𝑐2superscript𝜉23\displaystyle=\frac{2|\xi|^{2}\left(\omega^{2}+i\omega a-c^{2}|\xi|^{2}\right)% +8c^{2}|\xi|^{4}}{\left(\omega^{2}+i\omega a-c^{2}|\xi|^{2}\right)^{3}}.= divide start_ARG 2 | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) + 8 italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT end_ARG start_ARG ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT end_ARG .

Hence |2q2/c2|C/(1+|ξ|2)superscript2subscript𝑞2superscript𝑐2𝐶1superscript𝜉2|\partial^{2}q_{-2}/\partial c^{2}|\leq C/(1+|\xi|^{2})| ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT / ∂ italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | ≤ italic_C / ( 1 + | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) for some constant C>0𝐶0C>0italic_C > 0 that depends on ω𝜔\omegaitalic_ω, a𝑎aitalic_a and c𝑐citalic_c, but not on ξ𝜉\xiitalic_ξ. Hence, (28) is satisfied for n=2𝑛2n=-2italic_n = - 2 and Theorem 2 applies to our preconditioner.

We also have that

ΩQMP=ΩQP+Ω(QMQ)P=I+RΩ+Ω(QMQ)PsubscriptΩsubscript𝑄𝑀𝑃subscriptΩ𝑄𝑃subscriptΩsubscript𝑄𝑀𝑄𝑃𝐼subscript𝑅ΩsubscriptΩsubscript𝑄𝑀𝑄𝑃\displaystyle\mathcal{E}_{\Omega}Q_{M}P=\mathcal{E}_{\Omega}QP+\mathcal{E}_{% \Omega}(Q_{M}-Q)P=I+R_{\Omega}+\mathcal{E}_{\Omega}(Q_{M}-Q)Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P = caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P + caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT - italic_Q ) italic_P = italic_I + italic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT + caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT - italic_Q ) italic_P (31)

where the error operator RΩ=ΩQPI:Hs(Ω)Hs(Ω):subscript𝑅ΩsubscriptΩ𝑄𝑃𝐼superscript𝐻𝑠Ωsuperscript𝐻𝑠ΩR_{\Omega}=\mathcal{E}_{\Omega}QP-I:H^{s}(\Omega)\to H^{s}(\Omega)italic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT = caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P - italic_I : italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) is compact and the restriction operator Ω:Hs(d)Hs(Ω):subscriptΩsuperscript𝐻𝑠superscript𝑑superscript𝐻𝑠Ω\mathcal{E}_{\Omega}:H^{s}(\mathbb{R}^{d})\to H^{s}(\Omega)caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT : italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) is bounded. Both of these operators are described above Theorem 1. According to Theorem 2, we also have that

Ω(QMQ)PvHs(Ω)CM2PvHs2(Ω)CM2vHs(Ω)subscriptnormsubscriptΩsubscript𝑄𝑀𝑄𝑃𝑣superscript𝐻𝑠Ω𝐶superscript𝑀2subscriptnorm𝑃𝑣superscript𝐻𝑠2Ω𝐶superscript𝑀2subscriptnorm𝑣superscript𝐻𝑠Ω\displaystyle\|\mathcal{E}_{\Omega}(Q_{M}-Q)Pv\|_{H^{s}(\Omega)}\leq\frac{C}{M% ^{2}}\|Pv\|_{H^{s-2}(\Omega)}\leq\frac{C}{M^{2}}\|v\|_{H^{s}(\Omega)}∥ caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT - italic_Q ) italic_P italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) end_POSTSUBSCRIPT ≤ divide start_ARG italic_C end_ARG start_ARG italic_M start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ∥ italic_P italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s - 2 end_POSTSUPERSCRIPT ( roman_Ω ) end_POSTSUBSCRIPT ≤ divide start_ARG italic_C end_ARG start_ARG italic_M start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ∥ italic_v ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) end_POSTSUBSCRIPT (32)

for some constant C=C(Ω,c,a,ω)𝐶𝐶Ω𝑐𝑎𝜔C=C(\Omega,c,a,\omega)italic_C = italic_C ( roman_Ω , italic_c , italic_a , italic_ω ) independent of M𝑀Mitalic_M. Hence, the norm of Ω(QMQ)PsubscriptΩsubscript𝑄𝑀𝑄𝑃\mathcal{E}_{\Omega}(Q_{M}-Q)Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT ( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT - italic_Q ) italic_P can be made arbitrarily small by increasing the number of interpolation points M𝑀M\to\inftyitalic_M → ∞. Therefore, the interpolated preconditioned operator ΩQMPsubscriptΩsubscript𝑄𝑀𝑃\mathcal{E}_{\Omega}Q_{M}Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P is a small bounded perturbation of the Fredholm operator I+RΩ𝐼subscript𝑅ΩI+R_{\Omega}italic_I + italic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT. With these preliminaries, we are able to prove the following theorem regarding the ability of the interpolated preconditioner QMsubscript𝑄𝑀Q_{M}italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT to cluster the spectrum of the operator ΩQMPsubscriptΩsubscript𝑄𝑀𝑃\mathcal{E}_{\Omega}Q_{M}Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P.

Theorem 3.

Let Bρsubscript𝐵𝜌B_{\rho}italic_B start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT denote a disk in the complex plane centered at z=1𝑧1z=1italic_z = 1 with arbitrarily small radius ρ>0𝜌0\rho>0italic_ρ > 0. There exists M=M(ρ)𝑀𝑀𝜌M=M(\rho)italic_M = italic_M ( italic_ρ ) large enough such that the interpolated preconditioned operator ΩQMP:Hs(Ω)Hs(Ω):subscriptΩsubscript𝑄𝑀𝑃superscript𝐻𝑠Ωsuperscript𝐻𝑠Ω\mathcal{E}_{\Omega}Q_{M}P:H^{s}(\Omega)\to H^{s}(\Omega)caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P : italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) is Fredholm (with the same index as ΩQPsubscriptΩ𝑄𝑃\mathcal{E}_{\Omega}QPcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P), and its spectrum is contained in the disk Bρsubscript𝐵𝜌B_{\rho}italic_B start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT except for a finite number of its eigenvalues.

Proof.

This proof is a direct application of Theorem 3.16 in Chapter 4 of [30]. First, we note that, as a consequence of Theorem 1, the disk Bρsubscript𝐵𝜌B_{\rho}italic_B start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT contains all but a finite number of the eigenvalues of the Fredholm operator ΩQP=I+RΩsubscriptΩ𝑄𝑃𝐼subscript𝑅Ω\mathcal{E}_{\Omega}QP=I+R_{\Omega}caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P = italic_I + italic_R start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT. From (31)-(32), we have that

ΩQMPΩQPHs(Ω)Hs(Ω)CM2.subscriptnormsubscriptΩsubscript𝑄𝑀𝑃subscriptΩ𝑄𝑃superscript𝐻𝑠Ωsuperscript𝐻𝑠Ω𝐶superscript𝑀2\displaystyle\|\mathcal{E}_{\Omega}Q_{M}P-\mathcal{E}_{\Omega}QP\|_{H^{s}(% \Omega)\to H^{s}(\Omega)}\leq\frac{C}{M^{2}}.∥ caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P - caligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P ∥ start_POSTSUBSCRIPT italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) → italic_H start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT ( roman_Ω ) end_POSTSUBSCRIPT ≤ divide start_ARG italic_C end_ARG start_ARG italic_M start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG .

Hence, for M𝑀Mitalic_M sufficiently large, we can apply [30, Thm 5.17, Ch 4] to deduce that ΩQMPsubscriptΩsubscript𝑄𝑀𝑃\mathcal{E}_{\Omega}Q_{M}Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P is also Fredholm (with the dame index as ΩQPsubscriptΩ𝑄𝑃\mathcal{E}_{\Omega}QPcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q italic_P). Also, we apply [30, Thm 3.16, Ch 4] to assert that the disk Bρsubscript𝐵𝜌B_{\rho}italic_B start_POSTSUBSCRIPT italic_ρ end_POSTSUBSCRIPT contains the spectrum of ΩQMPsubscriptΩsubscript𝑄𝑀𝑃\mathcal{E}_{\Omega}Q_{M}Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P except for finitely many of its eigenvalues. ∎

The practical significance of Theorem 3 is similar to that of Theorem 1, namely, that the interpolated preconditioner QMsubscript𝑄𝑀Q_{M}italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT possesses the ability to cluster the spectrum of ΩQMPsubscriptΩsubscript𝑄𝑀𝑃\mathcal{E}_{\Omega}Q_{M}Pcaligraphic_E start_POSTSUBSCRIPT roman_Ω end_POSTSUBSCRIPT italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_P in the vicinity of z=1𝑧1z=1italic_z = 1 in the complex plane provided that the number of interpolation points M𝑀Mitalic_M is large enough. This property is needed for a preconditioner to accelerate the convergence of the GMRES iterative method as stated in [42, Part VI].

4.4 Absorbing layer

So far we have analyzed the proposed interpolated preconditioner QMsubscript𝑄𝑀Q_{M}italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT in the continuous setting where the Helmholtz equation is posed in all of dsuperscript𝑑\mathbb{R}^{d}blackboard_R start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT. However, in computational practice, this unbounded domain must be truncated to a bounded domain ΩΩ\Omegaroman_Ω and the effect of the Sommerfeld radiation condition should be incorporated into an absorbing layer/boundary condition. Here we offer a way to incorporate an absorbing layer into the definition of the preconditioner QMsubscript𝑄𝑀Q_{M}italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT. A common choice for achieving absorption is the complexification of the wave speed in the vicinity of the boundary of the truncated domain ΩΩ\Omegaroman_Ω [17, 38, 1]. Hence, the wave speed is replaced as follows

γ(x)=c(x)(1iζ(x))𝛾𝑥𝑐𝑥1𝑖𝜁𝑥\displaystyle\gamma(x)=c(x)(1-i\zeta(x))italic_γ ( italic_x ) = italic_c ( italic_x ) ( 1 - italic_i italic_ζ ( italic_x ) ) (33)

where c(x)𝑐𝑥c(x)italic_c ( italic_x ) is the real-valued wave speed, and ζ(x)𝜁𝑥\zeta(x)italic_ζ ( italic_x ) is a real-valued, non-negative, bounded function whose support defines the absorbing layer denoted by ΩALsubscriptΩAL\Omega_{\rm AL}roman_Ω start_POSTSUBSCRIPT roman_AL end_POSTSUBSCRIPT. Typically, 0ζ(x)10𝜁𝑥much-less-than10\leq\zeta(x)\ll 10 ≤ italic_ζ ( italic_x ) ≪ 1.

Recall from Subsection 4.1, that because c𝑐citalic_c is a real-valued function, we are able to interpolate the preconditioner Q𝑄Qitalic_Q over a one-dimensional interval [cmin,cmax]subscript𝑐minsubscript𝑐max[c_{\rm min},c_{\rm max}][ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ]. We wish to preserve this univariate interpolation property which is highly desirable from the computational point of view. Therefore, we assume that the wave speed is constant c(x)=co𝑐𝑥subscript𝑐𝑜c(x)=c_{o}italic_c ( italic_x ) = italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT in the absorbing layer ΩALsubscriptΩAL\Omega_{\rm AL}roman_Ω start_POSTSUBSCRIPT roman_AL end_POSTSUBSCRIPT. As a result, the complexified wave speed c(x)(1iζ(x))𝑐𝑥1𝑖𝜁𝑥c(x)(1-i\zeta(x))italic_c ( italic_x ) ( 1 - italic_i italic_ζ ( italic_x ) ) has values in the complex plane either in the horizontal line from cminsubscript𝑐minc_{\rm min}italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT to cmaxsubscript𝑐maxc_{\rm max}italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, or in the vertical line from cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT to co(1+iζmax)subscript𝑐𝑜1𝑖subscript𝜁maxc_{o}(1+i\zeta_{\rm max})italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ( 1 + italic_i italic_ζ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ) where ζmaxsubscript𝜁max\zeta_{\rm max}italic_ζ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT is the maximum value of the function ζ𝜁\zetaitalic_ζ. As before, we let {c1,c2,,cM}[cmin,cmax]subscript𝑐1subscript𝑐2subscript𝑐𝑀subscript𝑐minsubscript𝑐max\{c_{1},c_{2},...,c_{M}\}\subset[c_{\rm min},c_{\rm max}]{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } ⊂ [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] be a finite ordered set of distinct interpolation points such that c1=cminsubscript𝑐1subscript𝑐minc_{1}=c_{\rm min}italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT, cM=cmaxsubscript𝑐𝑀subscript𝑐maxc_{M}=c_{\rm max}italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT, and there is mosubscript𝑚𝑜m_{o}italic_m start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT such that cmo=cosubscript𝑐subscript𝑚𝑜subscript𝑐𝑜c_{m_{o}}=c_{o}italic_c start_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_POSTSUBSCRIPT = italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT. Also let {φ1,,φM}subscript𝜑1subscript𝜑𝑀\{\varphi_{1},...,\varphi_{M}\}{ italic_φ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_φ start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT } be an interpolating basis associated with the points {c1,c2,,cM}subscript𝑐1subscript𝑐2subscript𝑐𝑀\{c_{1},c_{2},...,c_{M}\}{ italic_c start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_c start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT }, that is, each φm:[cmin,cmax][0,1]:subscript𝜑𝑚subscript𝑐minsubscript𝑐max01\varphi_{m}:[c_{\rm min},c_{\rm max}]\to[0,1]italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT : [ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] → [ 0 , 1 ] is a function such that φm(cn)=0subscript𝜑𝑚subscript𝑐𝑛0\varphi_{m}(c_{n})=0italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) = 0 for all nm𝑛𝑚n\neq mitalic_n ≠ italic_m, and φm(cm)=1subscript𝜑𝑚subscript𝑐𝑚1\varphi_{m}(c_{m})=1italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) = 1 for each m=1,,M𝑚1𝑀m=1,...,Mitalic_m = 1 , … , italic_M. Similarly, let {ζ0,ζ1,ζ2,,ζM~}[0,ζmax]subscript𝜁0subscript𝜁1subscript𝜁2subscript𝜁~𝑀0subscript𝜁max\{\zeta_{0},\zeta_{1},\zeta_{2},...,\zeta_{\tilde{M}}\}\subset[0,\zeta_{\rm max}]{ italic_ζ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , italic_ζ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_ζ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , … , italic_ζ start_POSTSUBSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUBSCRIPT } ⊂ [ 0 , italic_ζ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] be a finite ordered set of distinct interpolation points such that ζ0=0subscript𝜁00\zeta_{0}=0italic_ζ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 and ζM~=ζmaxsubscript𝜁~𝑀subscript𝜁max\zeta_{\tilde{M}}=\zeta_{\rm max}italic_ζ start_POSTSUBSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUBSCRIPT = italic_ζ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT where M~=M~(M)~𝑀~𝑀𝑀\tilde{M}=\tilde{M}(M)over~ start_ARG italic_M end_ARG = over~ start_ARG italic_M end_ARG ( italic_M ). Also let {ϕ0,ϕ1,,ϕM~}subscriptitalic-ϕ0subscriptitalic-ϕ1subscriptitalic-ϕ~𝑀\{\phi_{0},\phi_{1},...,\phi_{\tilde{M}}\}{ italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , italic_ϕ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_ϕ start_POSTSUBSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUBSCRIPT } be an interpolating basis associated with the points {ζ0,ζ1,,ζM~}subscript𝜁0subscript𝜁1subscript𝜁~𝑀\{\zeta_{0},\zeta_{1},...,\zeta_{\tilde{M}}\}{ italic_ζ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT , italic_ζ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , … , italic_ζ start_POSTSUBSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUBSCRIPT }, that is, each ϕm:[coζmax,0][0,1]:subscriptitalic-ϕ𝑚subscript𝑐𝑜subscript𝜁max001\phi_{m}:[-c_{o}\zeta_{\rm max},0]\to[0,1]italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT : [ - italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT italic_ζ start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT , 0 ] → [ 0 , 1 ] is a function such that ϕm(coζn)=0subscriptitalic-ϕ𝑚subscript𝑐𝑜subscript𝜁𝑛0\phi_{m}(-c_{o}\zeta_{n})=0italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( - italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT italic_ζ start_POSTSUBSCRIPT italic_n end_POSTSUBSCRIPT ) = 0 for all nm𝑛𝑚n\neq mitalic_n ≠ italic_m, and ϕm(coζm)=1subscriptitalic-ϕ𝑚subscript𝑐𝑜subscript𝜁𝑚1\phi_{m}(-c_{o}\zeta_{m})=1italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( - italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT italic_ζ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) = 1 for each m=0,1,,M~𝑚01~𝑀m=0,1,...,\tilde{M}italic_m = 0 , 1 , … , over~ start_ARG italic_M end_ARG. Using these two interpolating bases, we define the pseudodifferential interpolating symbol as

qM(γ,ξ)=m=1Mφm(𝔢(γ))ϕ0(𝔪(γ))q(cm,ξ)+m=1M~φmo(𝔢(γ))ϕm(𝔪(γ))q(co(1iζm),ξ),subscript𝑞𝑀𝛾𝜉superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝔢𝛾subscriptitalic-ϕ0𝔪𝛾𝑞subscript𝑐𝑚𝜉superscriptsubscript𝑚1~𝑀subscript𝜑subscript𝑚𝑜𝔢𝛾subscriptitalic-ϕ𝑚𝔪𝛾𝑞subscript𝑐𝑜1𝑖subscript𝜁𝑚𝜉\displaystyle q_{M}(\gamma,\xi)=\sum_{m=1}^{M}\varphi_{m}(\mathfrak{Re}(\gamma% ))\phi_{0}(\mathfrak{Im}(\gamma))q(c_{m},\xi)+\sum_{m=1}^{\tilde{M}}\varphi_{m% _{o}}(\mathfrak{Re}(\gamma))\phi_{m}(\mathfrak{Im}(\gamma))q(c_{o}(1-i\zeta_{m% }),\xi),italic_q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT ( italic_γ , italic_ξ ) = ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( fraktur_R fraktur_e ( italic_γ ) ) italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( fraktur_I fraktur_m ( italic_γ ) ) italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) + ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( fraktur_R fraktur_e ( italic_γ ) ) italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( fraktur_I fraktur_m ( italic_γ ) ) italic_q ( italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ( 1 - italic_i italic_ζ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) , italic_ξ ) , (34)

and its associated pseudodifferential operator as

(QMv)(x)subscript𝑄𝑀𝑣𝑥\displaystyle(Q_{M}v)(x)( italic_Q start_POSTSUBSCRIPT italic_M end_POSTSUBSCRIPT italic_v ) ( italic_x ) =(2π)d/2m=1Mφm(𝔢(γ(x)))ϕ0(𝔪(γ(x)))q(cm,ξ)v^(ξ)eixξ𝑑ξabsentsuperscript2𝜋𝑑2superscriptsubscript𝑚1𝑀subscript𝜑𝑚𝔢𝛾𝑥subscriptitalic-ϕ0𝔪𝛾𝑥𝑞subscript𝑐𝑚𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle=(2\pi)^{-d/2}\sum_{m=1}^{M}\varphi_{m}(\mathfrak{Re}(\gamma(x)))% \phi_{0}(\mathfrak{Im}(\gamma(x)))\int q(c_{m},\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xi= ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_M end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( fraktur_R fraktur_e ( italic_γ ( italic_x ) ) ) italic_ϕ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( fraktur_I fraktur_m ( italic_γ ( italic_x ) ) ) ∫ italic_q ( italic_c start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ
+(2π)d/2m=1M~φmo(𝔢(γ(x)))ϕm(𝔪(γ(x)))q(co(1iζm),ξ)v^(ξ)eixξ𝑑ξ.superscript2𝜋𝑑2superscriptsubscript𝑚1~𝑀subscript𝜑subscript𝑚𝑜𝔢𝛾𝑥subscriptitalic-ϕ𝑚𝔪𝛾𝑥𝑞subscript𝑐𝑜1𝑖subscript𝜁𝑚𝜉^𝑣𝜉superscript𝑒𝑖𝑥𝜉differential-d𝜉\displaystyle\quad+(2\pi)^{-d/2}\sum_{m=1}^{\tilde{M}}\varphi_{m_{o}}(% \mathfrak{Re}(\gamma(x)))\phi_{m}(\mathfrak{Im}(\gamma(x)))\int q(c_{o}(1-i% \zeta_{m}),\xi)\hat{v}(\xi)e^{ix\cdot\xi}d\xi.+ ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∑ start_POSTSUBSCRIPT italic_m = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT over~ start_ARG italic_M end_ARG end_POSTSUPERSCRIPT italic_φ start_POSTSUBSCRIPT italic_m start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( fraktur_R fraktur_e ( italic_γ ( italic_x ) ) ) italic_ϕ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ( fraktur_I fraktur_m ( italic_γ ( italic_x ) ) ) ∫ italic_q ( italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ( 1 - italic_i italic_ζ start_POSTSUBSCRIPT italic_m end_POSTSUBSCRIPT ) , italic_ξ ) over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT italic_d italic_ξ . (35)

where q=q2𝑞subscript𝑞2q=q_{-2}italic_q = italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT is defined in (21).

5 Numerical experiments

In the numerical examples presented in this section, we consider a square domain ΩΩ\Omegaroman_Ω discretized with a regular grid with N𝑁Nitalic_N points in each coordinate direction. So there are N2superscript𝑁2N^{2}italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT evaluation points that represent the degrees of freedom. Let PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT denote the second-order centered finite difference discretization of the Helmholtz operator P𝑃Pitalic_P defined by (16) augmented by periodic boundary conditions. Also let QM,Nsubscript𝑄𝑀𝑁Q_{M,N}italic_Q start_POSTSUBSCRIPT italic_M , italic_N end_POSTSUBSCRIPT be the proposed pseudodifferential preconditioner (27) using the symbol q2subscript𝑞2q_{-2}italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT from (21), numerically evaluated using the FFT on the same N2superscript𝑁2N^{2}italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT-grid, and using M𝑀Mitalic_M interpolation terms as defined in Section 4. So the goal is to solve the following linear system

PNu=fNsubscript𝑃𝑁𝑢subscript𝑓𝑁\displaystyle P_{N}u=f_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT italic_u = italic_f start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT (36)

where PN:N×NN×N:subscript𝑃𝑁superscript𝑁𝑁superscript𝑁𝑁P_{N}:\mathbb{C}^{N\times N}\to\mathbb{C}^{N\times N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT : blackboard_C start_POSTSUPERSCRIPT italic_N × italic_N end_POSTSUPERSCRIPT → blackboard_C start_POSTSUPERSCRIPT italic_N × italic_N end_POSTSUPERSCRIPT and fNsubscript𝑓𝑁f_{N}italic_f start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT denotes the evaluation of (41) at the points on the N2superscript𝑁2N^{2}italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT-grid. The preconditioned version of (36) then becomes

(QN,MPN)u=QN,MfNsubscript𝑄𝑁𝑀subscript𝑃𝑁𝑢subscript𝑄𝑁𝑀subscript𝑓𝑁\displaystyle(Q_{N,M}P_{N})u=Q_{N,M}f_{N}( italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ) italic_u = italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_f start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT (37)

where the preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT is described above.

Our objective is to observe the computational complexity of the proposed preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT, its ability to precondition a discrete Helmholtz operator PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT, and observe the acceleration provided for the convergence of a GMRES solver for (37). The preconditioning performance will be characterized by the following spectral condition number,

cond(A)=maxλσ(A)|λ|minλσ(A)|λ|cond𝐴subscript𝜆𝜎𝐴𝜆subscript𝜆𝜎𝐴𝜆\displaystyle\text{cond}(A)=\frac{\max_{\lambda\in\sigma(A)}|\lambda|}{\min_{% \lambda\in\sigma(A)}|\lambda|}cond ( italic_A ) = divide start_ARG roman_max start_POSTSUBSCRIPT italic_λ ∈ italic_σ ( italic_A ) end_POSTSUBSCRIPT | italic_λ | end_ARG start_ARG roman_min start_POSTSUBSCRIPT italic_λ ∈ italic_σ ( italic_A ) end_POSTSUBSCRIPT | italic_λ | end_ARG (38)

where σ(A)𝜎𝐴\sigma(A)italic_σ ( italic_A ) denotes the spectrum of a generic matrix A𝐴Aitalic_A. Note that the spectral condition number (38) characterizes the standard condition number for normal (unitarily diagonalizable) matrices such as PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT obtained from the finite difference discretization of the Helmholtz operator. In all numerical experiments, the largest and smallest (in absolute value) eigenvalues were computed using MATLAB “eigs” function for a matrix-free calculations using function handles.

For comparison purposes, we also define a pseudodifferential preconditioner Q1subscript𝑄1Q_{1}italic_Q start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT with a single interpolation point at the background wave speed cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT as follows,

(Q1v)(x)=(2π)d/2v^(ξ)eixξω2+iωaoco2|ξ|2𝑑ξ,subscript𝑄1𝑣𝑥superscript2𝜋𝑑2^𝑣𝜉superscript𝑒𝑖𝑥𝜉superscript𝜔2𝑖𝜔subscript𝑎𝑜superscriptsubscript𝑐𝑜2superscript𝜉2differential-d𝜉\displaystyle(Q_{1}v)(x)=(2\pi)^{-d/2}\int\frac{\hat{v}(\xi)e^{ix\cdot\xi}}{% \omega^{2}+i\omega a_{o}-c_{o}^{2}|\xi|^{2}}d\xi,( italic_Q start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT italic_v ) ( italic_x ) = ( 2 italic_π ) start_POSTSUPERSCRIPT - italic_d / 2 end_POSTSUPERSCRIPT ∫ divide start_ARG over^ start_ARG italic_v end_ARG ( italic_ξ ) italic_e start_POSTSUPERSCRIPT italic_i italic_x ⋅ italic_ξ end_POSTSUPERSCRIPT end_ARG start_ARG italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_i italic_ω italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT - italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT | italic_ξ | start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG italic_d italic_ξ , (39)

along with its discrete counterpart QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT using the FFT for its numerical evaluation.

5.1 Example 1: Circular inclusion

We consider a square domain Ω2Ωsuperscript2\Omega\subset\mathbb{R}^{2}roman_Ω ⊂ blackboard_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT of side L=1𝐿1L=1italic_L = 1 centered at the origin. Periodic boundary conditions are assumed. The background medium has a wave speed co=1subscript𝑐𝑜1c_{o}=1italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT = 1. We place a circular inclusion of radius 0.050.050.050.05 with wavespeed c=co+δ𝑐subscript𝑐𝑜𝛿c=c_{o}+\deltaitalic_c = italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT + italic_δ. However, smooth transition is accomplished by the following definition,

c(x,y)=co+δHη(0.05x2+y2),𝑐𝑥𝑦subscript𝑐𝑜𝛿subscript𝐻𝜂0.05superscript𝑥2superscript𝑦2\displaystyle c(x,y)=c_{o}+\delta H_{\eta}(0.05-\sqrt{x^{2}+y^{2}}),italic_c ( italic_x , italic_y ) = italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT + italic_δ italic_H start_POSTSUBSCRIPT italic_η end_POSTSUBSCRIPT ( 0.05 - square-root start_ARG italic_x start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_y start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ) , (40)

where Hη(s)=1/(1+es/η)subscript𝐻𝜂𝑠11superscript𝑒𝑠𝜂H_{\eta}(s)=1/(1+e^{-s/\eta})italic_H start_POSTSUBSCRIPT italic_η end_POSTSUBSCRIPT ( italic_s ) = 1 / ( 1 + italic_e start_POSTSUPERSCRIPT - italic_s / italic_η end_POSTSUPERSCRIPT ) is a smooth version of the Heaviside function characterized by the transition parameter η>0𝜂0\eta>0italic_η > 0. The smaller η𝜂\etaitalic_η, the sharper the transition from 00 to 1111 in H(s)𝐻𝑠H(s)italic_H ( italic_s ), and consequently the sharper the transition from the background wave speed cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT to the inclusion wave speed co+δsubscript𝑐𝑜𝛿c_{o}+\deltaitalic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT + italic_δ. The source function f𝑓fitalic_f is defined as a plane wave

f(x,y)=eiω/cox𝑓𝑥𝑦superscript𝑒𝑖𝜔subscript𝑐𝑜𝑥\displaystyle f(x,y)=e^{i\omega/c_{o}x}italic_f ( italic_x , italic_y ) = italic_e start_POSTSUPERSCRIPT italic_i italic_ω / italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT italic_x end_POSTSUPERSCRIPT (41)

with wavenumber ω/co𝜔subscript𝑐𝑜\omega/c_{o}italic_ω / italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT where ω𝜔\omegaitalic_ω is the angular frequency of temporal oscillation. This experimental setup allows us to numerically test the proposed preconditioner under reproducible conditions and study its performance for a range of values for the number of interpolation points M𝑀Mitalic_M, the frequency ω𝜔\omegaitalic_ω, points per wavelength (PPW), the wave speed contrast parameter δ𝛿\deltaitalic_δ, the wave speed smoothness parameter η𝜂\etaitalic_η, and the damping coefficient a𝑎aitalic_a. The points per wavelength are defined as PPW =2πcoN/(Lω)absent2𝜋subscript𝑐𝑜𝑁𝐿𝜔=2\pi c_{o}N/(L\omega)= 2 italic_π italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT italic_N / ( italic_L italic_ω ).

Complexity:

First, we observe the complexity of evaluating the operator QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT. There results are displayed in Table 1 for N=180,360,720,1440,2880𝑁18036072014402880N=180,360,720,1440,2880italic_N = 180 , 360 , 720 , 1440 , 2880 and M=2,4,8𝑀248M=2,4,8italic_M = 2 , 4 , 8. Since the algorithm amounts to one FFT of the input, followed by M𝑀Mitalic_M inverse FFTs, then the assumed complexity in terms of CPU time is TNMNdβlogNsimilar-tosubscript𝑇𝑁𝑀superscript𝑁𝑑𝛽𝑁T_{N}\sim MN^{d\beta}\log Nitalic_T start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ∼ italic_M italic_N start_POSTSUPERSCRIPT italic_d italic_β end_POSTSUPERSCRIPT roman_log italic_N where d=2𝑑2d=2italic_d = 2 for a two-dimensional setting. Our objective is to empirically observe the order β𝛽\betaitalic_β to numerically test the estimate obtained in Lemma 1 which dictates that β=1𝛽1\beta=1italic_β = 1. An observed order of β=1𝛽1\beta=1italic_β = 1 would imply that, in practice, the proposed preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT possesses nearly linear complexity with respect to the DOF=N2superscript𝑁2N^{2}italic_N start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT in a two-dimensional setting. Table 1 shows that the order β𝛽\betaitalic_β is close to 1111 for each choice of the number of interpolation points M𝑀Mitalic_M in conformity with Lemma 1.

Table 1: Quantification of computational complexity through normalized CPU time for a range of grid points N𝑁Nitalic_N and interpolation nodes M𝑀Mitalic_M. The model for complexity is TNMNdβlogNsimilar-tosubscript𝑇𝑁𝑀superscript𝑁𝑑𝛽𝑁T_{N}\sim MN^{d\beta}\log Nitalic_T start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT ∼ italic_M italic_N start_POSTSUPERSCRIPT italic_d italic_β end_POSTSUPERSCRIPT roman_log italic_N where d=2𝑑2d=2italic_d = 2 for the two-dimensional setting. For each fixed M𝑀Mitalic_M, the observed order β𝛽\betaitalic_β is computed as the slope of a linear regression fitted through the CPU times TNsubscript𝑇𝑁T_{N}italic_T start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT versus N𝑁Nitalic_N in the log-log scale. A similar regression for all the values on the table renders an observed order β=1.03𝛽1.03\beta=1.03italic_β = 1.03. The CPU times displayed here are normalized to the CPU time for N=180𝑁180N=180italic_N = 180 and M=2𝑀2M=2italic_M = 2.
M𝑀Mitalic_M Normalized CPU Time Order β𝛽\betaitalic_β
N=180𝑁180N=180italic_N = 180 N=360𝑁360N=360italic_N = 360 N=720𝑁720N=720italic_N = 720 N=1440𝑁1440N=1440italic_N = 1440 N=2880𝑁2880N=2880italic_N = 2880
2222 1.00×1001.00superscript1001.00\times 10^{0}1.00 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 3.23×1003.23superscript1003.23\times 10^{0}3.23 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 1.43×1011.43superscript1011.43\times 10^{1}1.43 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.05×1016.05superscript1016.05\times 10^{1}6.05 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.62×1022.62superscript1022.62\times 10^{2}2.62 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.021.021.021.02
4444 1.62×1001.62superscript1001.62\times 10^{0}1.62 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 6.82×1006.82superscript1006.82\times 10^{0}6.82 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.67×1012.67superscript1012.67\times 10^{1}2.67 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.17×1021.17superscript1021.17\times 10^{2}1.17 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 5.06×1025.06superscript1025.06\times 10^{2}5.06 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.031.031.031.03
8888 3.18×1003.18superscript1003.18\times 10^{0}3.18 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 1.23×1011.23superscript1011.23\times 10^{1}1.23 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 5.45×1015.45superscript1015.45\times 10^{1}5.45 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.28×1022.28superscript1022.28\times 10^{2}2.28 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 9.98×1029.98superscript1029.98\times 10^{2}9.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.041.041.041.04

Frequency:

Here we explore the spectral behavior of the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT, for a range of frequencies ω𝜔\omegaitalic_ω and interpolation points M𝑀Mitalic_M, while holding the smoothness parameter η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, wave speed contrast parameter δ=4𝛿4\delta=4italic_δ = 4, damping coefficient a=20𝑎20a=20italic_a = 20, and points per wavelength PPW = 12121212 fixed. Table 2 displays the results from these experiments. We observe that PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT (finite difference discretization of the Helmholtz operator) is very ill-conditioned (cond(PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT) 10121017similar-toabsentsuperscript1012superscript1017\sim 10^{12}-10^{17}∼ 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT 17 end_POSTSUPERSCRIPT) and its condition number grows with the frequency ω𝜔\omegaitalic_ω. The preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point for the background wave speed cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT is able to reduce the condition number by several order of magnitude (cond(QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT) 107108similar-toabsentsuperscript107superscript108\sim 10^{7}-10^{8}∼ 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT) and behaves relatively stably with the frequency ω𝜔\omegaitalic_ω. The condition number of QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for M2𝑀2M\geq 2italic_M ≥ 2 is again several orders of magnitude smaller (cond(QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT) 101102similar-toabsentsuperscript101superscript102\sim 10^{1}-10^{2}∼ 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT). Additionally, we observe that the condition number of QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT decreases as the frequency ω𝜔\omegaitalic_ω increases. This behavior conforms with Remark 1 stating that the preconditioner Q𝑄Qitalic_Q inverts the Helmholtz operator P𝑃Pitalic_P up to an error of pseudodifferential order 11-1- 1. This means that the preconditioner becomes more effective for problems whose solutions are highly oscillatory (ω𝜔\omega\to\inftyitalic_ω → ∞) because the error operator behaves like a low-pass filter.

Table 2: Spectral condition number (38) for the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for increasing frequencies ω𝜔\omegaitalic_ω and interpolation points M𝑀Mitalic_M. For comparison, the preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point is also tested. In all cases, the wave speed smoothness η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, wave speed contrast parameter δ=4𝛿4\delta=4italic_δ = 4, damping coefficient a=20𝑎20a=20italic_a = 20, and points per wavelength PPW = 12121212 are fixed.
ω=20π𝜔20𝜋\omega=20\piitalic_ω = 20 italic_π ω=40π𝜔40𝜋\omega=40\piitalic_ω = 40 italic_π ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π ω=160π𝜔160𝜋\omega=160\piitalic_ω = 160 italic_π ω=320π𝜔320𝜋\omega=320\piitalic_ω = 320 italic_π
PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 7.95×10127.95superscript10127.95\times 10^{12}7.95 × 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT 1.31×10141.31superscript10141.31\times 10^{14}1.31 × 10 start_POSTSUPERSCRIPT 14 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 3.39×10163.39superscript10163.39\times 10^{16}3.39 × 10 start_POSTSUPERSCRIPT 16 end_POSTSUPERSCRIPT 5.42×10175.42superscript10175.42\times 10^{17}5.42 × 10 start_POSTSUPERSCRIPT 17 end_POSTSUPERSCRIPT
QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 1.21×1071.21superscript1071.21\times 10^{7}1.21 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 2.43×1072.43superscript1072.43\times 10^{7}2.43 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.69×1075.69superscript1075.69\times 10^{7}5.69 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 1.23×1081.23superscript1081.23\times 10^{8}1.23 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT 2.37×1082.37superscript1082.37\times 10^{8}2.37 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT
QN,2PNsubscript𝑄𝑁2subscript𝑃𝑁Q_{N,2}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 2 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 7.28×1027.28superscript1027.28\times 10^{2}7.28 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 4.23×1024.23superscript1024.23\times 10^{2}4.23 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.98×1021.98superscript1021.98\times 10^{2}1.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 7.79×1017.79superscript1017.79\times 10^{1}7.79 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 4.79×1014.79superscript1014.79\times 10^{1}4.79 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,4PNsubscript𝑄𝑁4subscript𝑃𝑁Q_{N,4}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 4 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 7.00×1027.00superscript1027.00\times 10^{2}7.00 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 2.07×1022.07superscript1022.07\times 10^{2}2.07 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.01×1021.01superscript1021.01\times 10^{2}1.01 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 4.32×1014.32superscript1014.32\times 10^{1}4.32 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.16×1012.16superscript1012.16\times 10^{1}2.16 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,8PNsubscript𝑄𝑁8subscript𝑃𝑁Q_{N,8}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 8 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 5.25×1025.25superscript1025.25\times 10^{2}5.25 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.83×1021.83superscript1021.83\times 10^{2}1.83 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 7.84×1017.84superscript1017.84\times 10^{1}7.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 3.27×1013.27superscript1013.27\times 10^{1}3.27 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.65×1011.65superscript1011.65\times 10^{1}1.65 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,16PNsubscript𝑄𝑁16subscript𝑃𝑁Q_{N,16}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 16 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.91×1024.91superscript1024.91\times 10^{2}4.91 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.40×1021.40superscript1021.40\times 10^{2}1.40 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 7.11×1017.11superscript1017.11\times 10^{1}7.11 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 3.05×1013.05superscript1013.05\times 10^{1}3.05 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.51×1011.51superscript1011.51\times 10^{1}1.51 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,32PNsubscript𝑄𝑁32subscript𝑃𝑁Q_{N,32}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 32 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.91×1024.91superscript1024.91\times 10^{2}4.91 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.32×1021.32superscript1021.32\times 10^{2}1.32 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 6.97×1016.97superscript1016.97\times 10^{1}6.97 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 3.17×1013.17superscript1013.17\times 10^{1}3.17 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.35×1011.35superscript1011.35\times 10^{1}1.35 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT

Points per wavelength:

Now we explore the spectral behavior of the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT, for a range of points per wavelength PPW and interpolation points M𝑀Mitalic_M, for fixed smoothness parameter η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, wave speed contrast parameter δ=4𝛿4\delta=4italic_δ = 4, damping coefficient a=20𝑎20a=20italic_a = 20, and frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π. Table 3 displays the results of these experiments. We observe a rapid increase in the condition number for the unconditioned matrix PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT which is tamed by the preconditioners QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT. The reduction in the condition number when using a single interpolation point (at the background wave speed cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT) is considerable, but using 2222 or more interpolation points brings down the condition number by several orders of magnitude more until a floor is reached. We note the robustness of the preconditioner, implying that when QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT is applied to PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT, the condition number remains practically constant as the PPW increase.

Table 3: Spectral condition number (38) for the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for points per wavelength PPW and interpolation points M𝑀Mitalic_M. For comparison, the preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point is also tested. In all cases, the wave speed smoothness η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, wave speed contrast parameter δ=4𝛿4\delta=4italic_δ = 4, damping coefficient a=20𝑎20a=20italic_a = 20, and frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π are fixed.
PPW =4absent4=4= 4 PPW =6absent6=6= 6 PPW =12absent12=12= 12 PPW =24absent24=24= 24 PPW =48absent48=48= 48
PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 7.62×10127.62superscript10127.62\times 10^{12}7.62 × 10 start_POSTSUPERSCRIPT 12 end_POSTSUPERSCRIPT 1.30×10141.30superscript10141.30\times 10^{14}1.30 × 10 start_POSTSUPERSCRIPT 14 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 3.39×10163.39superscript10163.39\times 10^{16}3.39 × 10 start_POSTSUPERSCRIPT 16 end_POSTSUPERSCRIPT 5.43×10175.43superscript10175.43\times 10^{17}5.43 × 10 start_POSTSUPERSCRIPT 17 end_POSTSUPERSCRIPT
QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 3.74×1073.74superscript1073.74\times 10^{7}3.74 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.02×1075.02superscript1075.02\times 10^{7}5.02 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.69×1075.69superscript1075.69\times 10^{7}5.69 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.90×1075.90superscript1075.90\times 10^{7}5.90 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.96×1075.96superscript1075.96\times 10^{7}5.96 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT
QN,2PNsubscript𝑄𝑁2subscript𝑃𝑁Q_{N,2}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 2 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 6.91×1016.91superscript1016.91\times 10^{1}6.91 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.37×1021.37superscript1021.37\times 10^{2}1.37 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.98×1021.98superscript1021.98\times 10^{2}1.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 2.19×1022.19superscript1022.19\times 10^{2}2.19 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 2.25×1022.25superscript1022.25\times 10^{2}2.25 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT
QN,4PNsubscript𝑄𝑁4subscript𝑃𝑁Q_{N,4}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 4 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.74×1014.74superscript1014.74\times 10^{1}4.74 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 8.32×1018.32superscript1018.32\times 10^{1}8.32 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.01×1021.01superscript1021.01\times 10^{2}1.01 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.09×1021.09superscript1021.09\times 10^{2}1.09 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.12×1021.12superscript1021.12\times 10^{2}1.12 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT
QN,8PNsubscript𝑄𝑁8subscript𝑃𝑁Q_{N,8}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 8 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.50×1014.50superscript1014.50\times 10^{1}4.50 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.40×1016.40superscript1016.40\times 10^{1}6.40 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.84×1017.84superscript1017.84\times 10^{1}7.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 8.20×1018.20superscript1018.20\times 10^{1}8.20 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 8.35×1018.35superscript1018.35\times 10^{1}8.35 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,16PNsubscript𝑄𝑁16subscript𝑃𝑁Q_{N,16}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 16 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.12×1014.12superscript1014.12\times 10^{1}4.12 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 5.90×1015.90superscript1015.90\times 10^{1}5.90 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.11×1017.11superscript1017.11\times 10^{1}7.11 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.62×1017.62superscript1017.62\times 10^{1}7.62 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.79×1017.79superscript1017.79\times 10^{1}7.79 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,32PNsubscript𝑄𝑁32subscript𝑃𝑁Q_{N,32}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 32 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 4.31×1014.31superscript1014.31\times 10^{1}4.31 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.08×1016.08superscript1016.08\times 10^{1}6.08 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.97×1016.97superscript1016.97\times 10^{1}6.97 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.40×1017.40superscript1017.40\times 10^{1}7.40 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.57×1017.57superscript1017.57\times 10^{1}7.57 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT

Wave speed contrast:

Now we test the performance of the preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT under a range of values for the wave speed contrast parameter η𝜂\etaitalic_η and its interaction with the number of interpolation points M𝑀Mitalic_M. Recall that the background wave speed is co=1subscript𝑐𝑜1c_{o}=1italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT = 1 and the circular inclusion has wave speed c=co+δ𝑐subscript𝑐𝑜𝛿c=c_{o}+\deltaitalic_c = italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT + italic_δ as defined in (40) with a smooth transition characterized by a fixed parameter η=1/800𝜂1800\eta=1/800italic_η = 1 / 800. The frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, damping coefficient a=20𝑎20a=20italic_a = 20, and points per wavelength PPW =12absent12=12= 12 are fixed too. Table 4 displays the spectral condition number for QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for these experiments. First, we observe that the larger the contrast in the wave speed, the larger the condition number of the discrete Helmholtz operator PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT, and the more need there is for the preconditioner. We also note that, for each fixed contrast parameter δ𝛿\deltaitalic_δ, the preconditioner improves as the number of interpolation points M𝑀Mitalic_M increases but a floor on its performance is quickly reached.

Table 4: Spectral condition number (38) for the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for various numbers of interpolation points M𝑀Mitalic_M and wave speed contrast parameters δ𝛿\deltaitalic_δ. For comparison, the preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point is also tested. In all cases, the wave speed smoothness η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, damping coefficient a=20𝑎20a=20italic_a = 20, and points per wavelength PPW = 12121212 are fixed.
δ=1𝛿1\delta=1italic_δ = 1 δ=2𝛿2\delta=2italic_δ = 2 δ=4𝛿4\delta=4italic_δ = 4 δ=8𝛿8\delta=8italic_δ = 8 δ=16𝛿16\delta=16italic_δ = 16
PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 5.33×10135.33superscript10135.33\times 10^{13}5.33 × 10 start_POSTSUPERSCRIPT 13 end_POSTSUPERSCRIPT 2.72×10142.72superscript10142.72\times 10^{14}2.72 × 10 start_POSTSUPERSCRIPT 14 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.22×10162.22superscript10162.22\times 10^{16}2.22 × 10 start_POSTSUPERSCRIPT 16 end_POSTSUPERSCRIPT 2.83×10172.83superscript10172.83\times 10^{17}2.83 × 10 start_POSTSUPERSCRIPT 17 end_POSTSUPERSCRIPT
QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.70×1059.70superscript1059.70\times 10^{5}9.70 × 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT 6.47×1066.47superscript1066.47\times 10^{6}6.47 × 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT 5.69×1075.69superscript1075.69\times 10^{7}5.69 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 6.28×1086.28superscript1086.28\times 10^{8}6.28 × 10 start_POSTSUPERSCRIPT 8 end_POSTSUPERSCRIPT 8.14×1098.14superscript1098.14\times 10^{9}8.14 × 10 start_POSTSUPERSCRIPT 9 end_POSTSUPERSCRIPT
QN,2PNsubscript𝑄𝑁2subscript𝑃𝑁Q_{N,2}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 2 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.54×1009.54superscript1009.54\times 10^{0}9.54 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 3.22×1013.22superscript1013.22\times 10^{1}3.22 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.98×1021.98superscript1021.98\times 10^{2}1.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.90×1031.90superscript1031.90\times 10^{3}1.90 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT 2.37×1042.37superscript1042.37\times 10^{4}2.37 × 10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT
QN,4PNsubscript𝑄𝑁4subscript𝑃𝑁Q_{N,4}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 4 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 1.01×1011.01superscript1011.01\times 10^{1}1.01 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.60×1012.60superscript1012.60\times 10^{1}2.60 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.01×1021.01superscript1021.01\times 10^{2}1.01 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 5.59×1025.59superscript1025.59\times 10^{2}5.59 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 6.31×1036.31superscript1036.31\times 10^{3}6.31 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT
QN,8PNsubscript𝑄𝑁8subscript𝑃𝑁Q_{N,8}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 8 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.86×1009.86superscript1009.86\times 10^{0}9.86 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.35×1012.35superscript1012.35\times 10^{1}2.35 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.84×1017.84superscript1017.84\times 10^{1}7.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.71×1022.71superscript1022.71\times 10^{2}2.71 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.55×1031.55superscript1031.55\times 10^{3}1.55 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT
QN,16PNsubscript𝑄𝑁16subscript𝑃𝑁Q_{N,16}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 16 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.66×1009.66superscript1009.66\times 10^{0}9.66 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.34×1012.34superscript1012.34\times 10^{1}2.34 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.11×1017.11superscript1017.11\times 10^{1}7.11 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.36×1022.36superscript1022.36\times 10^{2}2.36 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.08×1031.08superscript1031.08\times 10^{3}1.08 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT
QN,32PNsubscript𝑄𝑁32subscript𝑃𝑁Q_{N,32}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 32 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.64×1009.64superscript1009.64\times 10^{0}9.64 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.33×1012.33superscript1012.33\times 10^{1}2.33 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.97×1016.97superscript1016.97\times 10^{1}6.97 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.59×1022.59superscript1022.59\times 10^{2}2.59 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.06×1031.06superscript1031.06\times 10^{3}1.06 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT

Wave speed smoothness:

Here we test the performance of the preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT for a range of values for the smoothness parameter η𝜂\etaitalic_η and its interaction with the number of interpolation points M𝑀Mitalic_M. The other parameters δ=4𝛿4\delta=4italic_δ = 4, ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, a=20𝑎20a=20italic_a = 20, and PPW =12absent12=12= 12 are fixed. The results from these experiments are shown in Table 5. As expected from the mathematical theory of pseudodifferential operators and the accuracy of interpolation, we observe that the proposed preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT is more effective when the wave speed profile is smoother. As before, when using a single interpolation point (at the background wave speed cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT), the reduction in the condition number is significant, but using 2222 or more interpolation points brings down the condition number by several orders of magnitude more.

Table 5: Spectral condition number (38) for the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for various numbers of interpolation points M𝑀Mitalic_M and wave speed smoothness parameters η𝜂\etaitalic_η. For comparison, the preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point is also tested. In all cases, the wave speed contrast δ=4𝛿4\delta=4italic_δ = 4, frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, damping coefficient a=20𝑎20a=20italic_a = 20, and points per wavelength PPW = 12121212 are fixed.
η=1/200𝜂1200\eta=1/200italic_η = 1 / 200 η=1/400𝜂1400\eta=1/400italic_η = 1 / 400 η=1/800𝜂1800\eta=1/800italic_η = 1 / 800 η=1/1600𝜂11600\eta=1/1600italic_η = 1 / 1600 η=1/3200𝜂13200\eta=1/3200italic_η = 1 / 3200
PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 2.10×10152.10superscript10152.10\times 10^{15}2.10 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT
QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 2.09×1072.09superscript1072.09\times 10^{7}2.09 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 4.11×1074.11superscript1074.11\times 10^{7}4.11 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.69×1075.69superscript1075.69\times 10^{7}5.69 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 6.44×1076.44superscript1076.44\times 10^{7}6.44 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 6.72×1076.72superscript1076.72\times 10^{7}6.72 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT
QN,2PNsubscript𝑄𝑁2subscript𝑃𝑁Q_{N,2}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 2 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 8.15×1018.15superscript1018.15\times 10^{1}8.15 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.87×1016.87superscript1016.87\times 10^{1}6.87 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.98×1021.98superscript1021.98\times 10^{2}1.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 6.87×1026.87superscript1026.87\times 10^{2}6.87 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.76×1031.76superscript1031.76\times 10^{3}1.76 × 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT
QN,4PNsubscript𝑄𝑁4subscript𝑃𝑁Q_{N,4}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 4 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 7.93×1007.93superscript1007.93\times 10^{0}7.93 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 3.08×1013.08superscript1013.08\times 10^{1}3.08 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.01×1021.01superscript1021.01\times 10^{2}1.01 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 3.03×1023.03superscript1023.03\times 10^{2}3.03 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 9.80×1029.80superscript1029.80\times 10^{2}9.80 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT
QN,8PNsubscript𝑄𝑁8subscript𝑃𝑁Q_{N,8}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 8 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.55×1009.55superscript1009.55\times 10^{0}9.55 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.43×1012.43superscript1012.43\times 10^{1}2.43 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.84×1017.84superscript1017.84\times 10^{1}7.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.13×1022.13superscript1022.13\times 10^{2}2.13 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 5.77×1025.77superscript1025.77\times 10^{2}5.77 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT
QN,16PNsubscript𝑄𝑁16subscript𝑃𝑁Q_{N,16}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 16 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.41×1009.41superscript1009.41\times 10^{0}9.41 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.36×1012.36superscript1012.36\times 10^{1}2.36 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.11×1017.11superscript1017.11\times 10^{1}7.11 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 2.24×1022.24superscript1022.24\times 10^{2}2.24 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 5.83×1025.83superscript1025.83\times 10^{2}5.83 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT
QN,32PNsubscript𝑄𝑁32subscript𝑃𝑁Q_{N,32}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 32 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 8.61×1008.61superscript1008.61\times 10^{0}8.61 × 10 start_POSTSUPERSCRIPT 0 end_POSTSUPERSCRIPT 2.50×1012.50superscript1012.50\times 10^{1}2.50 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.97×1016.97superscript1016.97\times 10^{1}6.97 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 1.92×1021.92superscript1021.92\times 10^{2}1.92 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 4.95×1024.95superscript1024.95\times 10^{2}4.95 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT

Damping:

We also test the performance of the preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT for a range of values for the damping coefficient a𝑎aitalic_a and number of interpolation points M𝑀Mitalic_M. The other parameters δ=4𝛿4\delta=4italic_δ = 4, η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, and PPW =12absent12=12= 12 are fixed. The results from these experiments are displayed in Table 6. We observe that the smaller the damping coefficient, the less effective the preconditioner QN,Msubscript𝑄𝑁𝑀Q_{N,M}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT. However, over this range of value for the damping coefficient, the preconditioner seems to be quite robust.

Table 6: Spectral condition number (38) for the discrete operators PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT and QN,MPNsubscript𝑄𝑁𝑀subscript𝑃𝑁Q_{N,M}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , italic_M end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT for various interpolation points M𝑀Mitalic_M and damping coefficients a𝑎aitalic_a. For comparison, the preconditioner QN,1subscript𝑄𝑁1Q_{N,1}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT with a single interpolation point is also tested. In all cases, the wave speed contrast δ=4𝛿4\delta=4italic_δ = 4, wave speed smoothness η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, frequency ω=80π𝜔80𝜋\omega=80\piitalic_ω = 80 italic_π, and points per wavelength PPW = 12121212 are fixed.
a=5𝑎5a=5italic_a = 5 a=10𝑎10a=10italic_a = 10 a=20𝑎20a=20italic_a = 20 a=40𝑎40a=40italic_a = 40 a=80𝑎80a=80italic_a = 80
PNsubscript𝑃𝑁P_{N}italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT 2.11×10152.11superscript10152.11\times 10^{15}2.11 × 10 start_POSTSUPERSCRIPT 15 end_POSTSUPERSCRIPT
QN,1PNsubscript𝑄𝑁1subscript𝑃𝑁Q_{N,1}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 1 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 8.92×1078.92superscript1078.92\times 10^{7}8.92 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 6.83×1076.83superscript1076.83\times 10^{7}6.83 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 5.69×1075.69superscript1075.69\times 10^{7}5.69 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 4.12×1074.12superscript1074.12\times 10^{7}4.12 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT 2.47×1072.47superscript1072.47\times 10^{7}2.47 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT
QN,2PNsubscript𝑄𝑁2subscript𝑃𝑁Q_{N,2}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 2 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 2.74×1022.74superscript1022.74\times 10^{2}2.74 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 2.42×1022.42superscript1022.42\times 10^{2}2.42 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.98×1021.98superscript1021.98\times 10^{2}1.98 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.41×1021.41superscript1021.41\times 10^{2}1.41 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 8.84×1018.84superscript1018.84\times 10^{1}8.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,4PNsubscript𝑄𝑁4subscript𝑃𝑁Q_{N,4}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 4 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 1.69×1021.69superscript1021.69\times 10^{2}1.69 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.25×1021.25superscript1021.25\times 10^{2}1.25 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 1.01×1021.01superscript1021.01\times 10^{2}1.01 × 10 start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT 8.09×1018.09superscript1018.09\times 10^{1}8.09 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 5.81×1015.81superscript1015.81\times 10^{1}5.81 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,8PNsubscript𝑄𝑁8subscript𝑃𝑁Q_{N,8}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 8 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 9.99×1019.99superscript1019.99\times 10^{1}9.99 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 8.40×1018.40superscript1018.40\times 10^{1}8.40 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.84×1017.84superscript1017.84\times 10^{1}7.84 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.61×1016.61superscript1016.61\times 10^{1}6.61 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 4.81×1014.81superscript1014.81\times 10^{1}4.81 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,16PNsubscript𝑄𝑁16subscript𝑃𝑁Q_{N,16}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 16 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 8.49×1018.49superscript1018.49\times 10^{1}8.49 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.81×1017.81superscript1017.81\times 10^{1}7.81 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.11×1017.11superscript1017.11\times 10^{1}7.11 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.07×1016.07superscript1016.07\times 10^{1}6.07 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 4.62×1014.62superscript1014.62\times 10^{1}4.62 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT
QN,32PNsubscript𝑄𝑁32subscript𝑃𝑁Q_{N,32}P_{N}italic_Q start_POSTSUBSCRIPT italic_N , 32 end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT italic_N end_POSTSUBSCRIPT 8.30×1018.30superscript1018.30\times 10^{1}8.30 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 7.54×1017.54superscript1017.54\times 10^{1}7.54 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.97×1016.97superscript1016.97\times 10^{1}6.97 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 6.08×1016.08superscript1016.08\times 10^{1}6.08 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT 4.61×1014.61superscript1014.61\times 10^{1}4.61 × 10 start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT

GMRES implementation:

Finally we also report some numerical results for solving the Helmholtz equation with and without the proposed preconditioner (27) using the GMRES iterative method. MATLAB “gmres” function was employed for matrix-free calculations using function handles. The restart parameter was set equal to 10101010. For the problem described by (40)-(41), the evolution of the residuals over the GMRES iterations are illustrated in Figure 1 for frequencies ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π, ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π and ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π, which respectively corresponds to 200200200200, 300300300300 and 400400400400 wavelengths across the domain. In all cases, the damping coefficient a=20𝑎20a=20italic_a = 20, the wave speed contrast δ=1𝛿1\delta=1italic_δ = 1, the wave speed smoothness η=1/800𝜂1800\eta=1/800italic_η = 1 / 800, and the number of interpolation points is M=8𝑀8M=8italic_M = 8 are fixed. The number of points per wavelength is PPW = 12121212 which leads to N=2400𝑁2400N=2400italic_N = 2400, N=3600𝑁3600N=3600italic_N = 3600 and N=4800𝑁4800N=4800italic_N = 4800, respectively. Hence, DOF = 5.760×1065.760superscript1065.760\times 10^{6}5.760 × 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT, DOF = 1.296×1071.296superscript1071.296\times 10^{7}1.296 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT and DOF = 2.304×1072.304superscript1072.304\times 10^{7}2.304 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT, respectively.

We observe from Figure 1 that applying the GMRES method naively to (36) leads to a very slow convergence pattern. By contrast, the application of the GMRES method to the preconditioned system (37) leads to much faster convergence. In fact, at 100100100100 iterations, the residual is at least 7777 orders of magnitude smaller when the proposed preconditioner is employed. For reference, the amplitude of the numerical solutions obtained after 100100100100 iterations are displayed in Figure 1.

Refer to caption
(a) Amplitude of the solution for ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π
Refer to caption
(b) GMRES residuals for ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π
Refer to caption
(c) Amplitude of the solution for ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π
Refer to caption
(d) GMRES residuals for ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π
Refer to caption
(e) Amplitude of the solution for ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π
Refer to caption
(f) GMRES residuals for ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π
Figure 1: Amplitude of the solution u𝑢uitalic_u (normalized in the Lsuperscript𝐿L^{\infty}italic_L start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT-norm) and relative residual for the first 100 iterations of the GMRES method applied to (36) and to (37) for frequencies ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π, ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π and ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π, respectively. In all cases, the GMRES+Predonditioner method reaches residuals more than 7777 orders of magnitude smaller than the GMRES method at the 100100100100th iteration. The wave speed profile c𝑐citalic_c and source f𝑓fitalic_f are defined by (40) (using δ=1𝛿1\delta=1italic_δ = 1 and η=1/800𝜂1800\eta=1/800italic_η = 1 / 800) and (41), respectively.

5.2 Example 2: Human head phantom and absorbing layer

In this subsection we apply the proposed preconditioner to solve the Helmholtz equation in a domain whose acoustic properties conform to those of a human head as defined by a synthetic phantom available as supplemental material from [6]. The spatial dimensions of the phantom are scaled in order to fit in a square domain Ω2Ωsuperscript2\Omega\subset\mathbb{R}^{2}roman_Ω ⊂ blackboard_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT of side L=1𝐿1L=1italic_L = 1 centered at the origin. Periodic boundary conditions are assumed at the sides of the square ΩΩ\Omegaroman_Ω, but an absorbing layer is introduced to enforce the outgoing behavior of the waves. This absorbing layer follows the setup described in Subsection 4.4. Specifically, the function ζ𝜁\zetaitalic_ζ defining the absorbing layer is given by

ζ(x1,x2)={0,if x12+x22<0.375,0.4(x12+x220.375),if x12+x220.375.𝜁subscript𝑥1subscript𝑥2cases0if x12+x22<0.375,0.4superscriptsubscript𝑥12superscriptsubscript𝑥220.375if x12+x220.375.\displaystyle\zeta(x_{1},x_{2})=\begin{cases}0,&\text{if $\sqrt{x_{1}^{2}+x_{2% }^{2}}<0.375$,}\\ 0.4(\sqrt{x_{1}^{2}+x_{2}^{2}}-0.375),&\text{if $\sqrt{x_{1}^{2}+x_{2}^{2}}% \geq 0.375$.}\end{cases}italic_ζ ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) = { start_ROW start_CELL 0 , end_CELL start_CELL if square-root start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG < 0.375 , end_CELL end_ROW start_ROW start_CELL 0.4 ( square-root start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG - 0.375 ) , end_CELL start_CELL if square-root start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ≥ 0.375 . end_CELL end_ROW (42)

The background medium has a wave speed co=1subscript𝑐𝑜1c_{o}=1italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT = 1 and damping coefficient ao=20.9subscript𝑎𝑜20.9a_{o}=20.9italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT = 20.9. The human head phantom also contains a cross section of the skull whose wave speed is cskull=6.71subscript𝑐skull6.71c_{\rm skull}=6.71italic_c start_POSTSUBSCRIPT roman_skull end_POSTSUBSCRIPT = 6.71 and damping coefficient askull=258subscript𝑎skull258a_{\rm skull}=258italic_a start_POSTSUBSCRIPT roman_skull end_POSTSUBSCRIPT = 258. The soft tissues representing the brain matter and fluids contained by the skull are assumed to have a uniform wave speed cbrain=4.55subscript𝑐brain4.55c_{\rm brain}=4.55italic_c start_POSTSUBSCRIPT roman_brain end_POSTSUBSCRIPT = 4.55 and damping coefficient abrain=10.7subscript𝑎brain10.7a_{\rm brain}=10.7italic_a start_POSTSUBSCRIPT roman_brain end_POSTSUBSCRIPT = 10.7.

The source function f𝑓fitalic_f is defined as a plane wave

f=(ω2(co2c2)co2(1+iao/ω)+iω(aao))uinc𝑓superscript𝜔2superscriptsubscript𝑐𝑜2superscript𝑐2superscriptsubscript𝑐𝑜21𝑖subscript𝑎𝑜𝜔𝑖𝜔𝑎subscript𝑎𝑜subscript𝑢inc\displaystyle f=-\left(\omega^{2}\frac{(c_{o}^{2}-c^{2})}{c_{o}^{2}}(1+ia_{o}/% \omega)+i\omega(a-a_{o})\right)u_{\rm inc}italic_f = - ( italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT divide start_ARG ( italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_c start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG ( 1 + italic_i italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT / italic_ω ) + italic_i italic_ω ( italic_a - italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ) ) italic_u start_POSTSUBSCRIPT roman_inc end_POSTSUBSCRIPT (43)

where ω𝜔\omegaitalic_ω is the angular frequency of temporal oscillation, and the incident wave field uincsubscript𝑢incu_{\rm inc}italic_u start_POSTSUBSCRIPT roman_inc end_POSTSUBSCRIPT is given by a Gaussian beam traveling downward (negative direction on the x2subscript𝑥2x_{2}italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT-axis),

uinc=WoW(x2x2,o)exp([x1x1,oW(x2x2,o)]2)exp(ik[x2x2,o+12κ(x2)(x1x1,o)2])subscript𝑢incsubscript𝑊𝑜𝑊subscript𝑥2subscript𝑥2𝑜superscriptdelimited-[]subscript𝑥1subscript𝑥1𝑜𝑊subscript𝑥2subscript𝑥2𝑜2𝑖𝑘delimited-[]subscript𝑥2subscript𝑥2𝑜12𝜅subscript𝑥2superscriptsubscript𝑥1subscript𝑥1𝑜2\displaystyle u_{\rm inc}=\frac{W_{o}}{W(x_{2}-x_{2,o})}\exp\left(-\left[\frac% {x_{1}-x_{1,o}}{W(x_{2}-x_{2,o})}\right]^{2}\right)\exp\left(-ik\left[x_{2}-x_% {2,o}+\frac{1}{2}\kappa(x_{2})(x_{1}-x_{1,o})^{2}\right]\right)italic_u start_POSTSUBSCRIPT roman_inc end_POSTSUBSCRIPT = divide start_ARG italic_W start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_ARG start_ARG italic_W ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT 2 , italic_o end_POSTSUBSCRIPT ) end_ARG roman_exp ( - [ divide start_ARG italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT 1 , italic_o end_POSTSUBSCRIPT end_ARG start_ARG italic_W ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT 2 , italic_o end_POSTSUBSCRIPT ) end_ARG ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) roman_exp ( - italic_i italic_k [ italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT 2 , italic_o end_POSTSUBSCRIPT + divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_κ ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) ( italic_x start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_x start_POSTSUBSCRIPT 1 , italic_o end_POSTSUBSCRIPT ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ] ) (44)
whereW(x2)=Wo1+(x2xR)2,k=ωco[1+iaoω]1/2,andκ(x2)=x2x22+xR2.formulae-sequencewhere𝑊subscript𝑥2subscript𝑊𝑜1superscriptsubscript𝑥2subscript𝑥R2formulae-sequence𝑘𝜔subscript𝑐𝑜superscriptdelimited-[]1𝑖subscript𝑎𝑜𝜔12and𝜅subscript𝑥2subscript𝑥2superscriptsubscript𝑥22superscriptsubscript𝑥R2\displaystyle\text{where}\quad W(x_{2})=W_{o}\sqrt{1+\left(\frac{x_{2}}{x_{\rm R% }}\right)^{2}},\quad k=\frac{\omega}{c_{o}}\left[1+i\frac{a_{o}}{\omega}\right% ]^{1/2},\quad\text{and}\quad\kappa(x_{2})=\frac{x_{2}}{x_{2}^{2}+x_{\rm R}^{2}}.where italic_W ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) = italic_W start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT square-root start_ARG 1 + ( divide start_ARG italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_x start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_k = divide start_ARG italic_ω end_ARG start_ARG italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_ARG [ 1 + italic_i divide start_ARG italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT end_ARG start_ARG italic_ω end_ARG ] start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT , and italic_κ ( italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) = divide start_ARG italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG italic_x start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_x start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG . (45)

Here, the point (x1,o,x2,o)=(0,0.3)subscript𝑥1𝑜subscript𝑥2𝑜00.3(x_{1,o},x_{2,o})=(0,0.3)( italic_x start_POSTSUBSCRIPT 1 , italic_o end_POSTSUBSCRIPT , italic_x start_POSTSUBSCRIPT 2 , italic_o end_POSTSUBSCRIPT ) = ( 0 , 0.3 ) is the focus of the beam and Wo=0.01subscript𝑊𝑜0.01W_{o}=0.01italic_W start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT = 0.01 is its waist, cosubscript𝑐𝑜c_{o}italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT is the background wave speed, aosubscript𝑎𝑜a_{o}italic_a start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT is the background damping coefficient, κ(y)𝜅𝑦\kappa(y)italic_κ ( italic_y ) is the curvature of the beam, and yR=ωWo2/(2co)subscript𝑦R𝜔superscriptsubscript𝑊𝑜22subscript𝑐𝑜y_{\rm R}=\omega W_{o}^{2}/(2c_{o})italic_y start_POSTSUBSCRIPT roman_R end_POSTSUBSCRIPT = italic_ω italic_W start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / ( 2 italic_c start_POSTSUBSCRIPT italic_o end_POSTSUBSCRIPT ) is called the Rayleigh range which is the point at which the curvature of the beam is the largest. The source f𝑓fitalic_f given by (43) in terms of the incident field uincsubscript𝑢incu_{\rm inc}italic_u start_POSTSUBSCRIPT roman_inc end_POSTSUBSCRIPT defined by (44)-(45) corresponds to a scattering problem for an acoustic beam propagating in the background medium and impinging on the model of the human head. Therefore, the solution to (1) for this source f𝑓fitalic_f corresponds to an acoustic field scattered by the human head acting as a penetrable obstacle. The (real-valued) wave speed c𝑐citalic_c and the function ζ𝜁\zetaitalic_ζ defining the absorbing layer are illustrated in Figure 2

Refer to caption
(a) Real-valued wave speed c𝑐citalic_c
Refer to caption
(b) Function ζ𝜁\zetaitalic_ζ for absorbing layer
Figure 2: (a) Wave speed profile c𝑐citalic_c from the supplemental material of [6]. (b) Function ζ𝜁\zetaitalic_ζ for the absorbing layer given by (42). The skull contour is shown for reference. The complex-valued wave speed is defined by (33).

For the problem described by (42)-(45), the evolution of the residuals over the GMRES iterations are illustrated in Figure 3 for frequencies ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π, ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π and ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π, which respectively corresponds to 200200200200, 300300300300 and 400400400400 wavelengths across the domain. In all cases, the the number of interpolation points is M=4𝑀4M=4italic_M = 4 for the real values of the wave speed and M~=3~𝑀3\tilde{M}=3over~ start_ARG italic_M end_ARG = 3 for the imaginary values of the wave speed (absorbing layer). The number of points per wavelength is PPW = 12121212 which leads to N=2400𝑁2400N=2400italic_N = 2400, N=3600𝑁3600N=3600italic_N = 3600 and N=4800𝑁4800N=4800italic_N = 4800, respectively. Hence, DOF = 5.760×1065.760superscript1065.760\times 10^{6}5.760 × 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT, DOF = 1.296×1071.296superscript1071.296\times 10^{7}1.296 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT and DOF = 2.304×1072.304superscript1072.304\times 10^{7}2.304 × 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT, respectively. As before, we used the MATLAB “gmres” function for matrix-free calculations using function handles. The restart parameter was set equal to 10101010.

We observe from Figure 3 that applying the GMRES method naively to (36) leads to a slow convergence pattern and that the application of the GMRES method to the preconditioned system (37) improves the convergence considerably. In fact, at 100100100100 iterations, the residual is at least 8888 orders of magnitude smaller when the proposed preconditioner is employed. For reference, the (log base 10) amplitude of the numerical solutions obtained after 100100100100 iterations are displayed in Figure 3.

Refer to caption
(a) Solution amplitude for ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π (log-scale)
Refer to caption
(b) Residuals for ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π
Refer to caption
(c) Solution amplitude for ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π (log-scale)
Refer to caption
(d) Residuals for ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π
Refer to caption
(e) Solution amplitude for ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π (log-scale)
Refer to caption
(f) Residuals for ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π
Figure 3: Amplitude of the solution u𝑢uitalic_u (log-scale) and relative residual for the first 100 iterations of the GMRES method applied to (36) and to (37) for frequencies ω=400π𝜔400𝜋\omega=400\piitalic_ω = 400 italic_π, ω=600π𝜔600𝜋\omega=600\piitalic_ω = 600 italic_π and ω=800π𝜔800𝜋\omega=800\piitalic_ω = 800 italic_π, respectively.


6 Conclusion

In this paper, we proposed a new pseudodifferential preconditioner for the Helmholtz equation in variable media with absorption. This preconditioner is constructed from the principal symbol q2=q2(x,ξ)subscript𝑞2subscript𝑞2𝑥𝜉q_{-2}=q_{-2}(x,\xi)italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT = italic_q start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT ( italic_x , italic_ξ ) (see (21)) of a non-classical pseudodifferential expansion for the inverse of the Helmholtz operator. We take advantage of the explicit dependence of this principal symbol on the wave speed c𝑐citalic_c to interpolate it as a function of the wave speed c𝑐citalic_c rather than as a function of the phase-space variables (x,ξ)𝑥𝜉(x,\xi)( italic_x , italic_ξ ). As a result we obtain an accurate and computationally efficient approximation of this symbol using univariate interpolation theory even when the original wave problem is posed in two- or three-dimensional domains. The overall complexity of the proposed preconditioner is 𝒪(MNdlogN)𝒪𝑀superscript𝑁𝑑𝑁\mathcal{O}(MN^{d}\log N)caligraphic_O ( italic_M italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT roman_log italic_N ) where M𝑀Mitalic_M is the number of interpolation points sampling the range of the wave speed, Ndsuperscript𝑁𝑑N^{d}italic_N start_POSTSUPERSCRIPT italic_d end_POSTSUPERSCRIPT is the DOF, and d=2,3𝑑23d=2,3italic_d = 2 , 3 is the dimension of the physical domain. The number M𝑀Mitalic_M of interpolation points is independent of the DOF, of the frequency ω𝜔\omegaitalic_ω, and of the dimension d𝑑ditalic_d. The proposed preconditioner can be implemented as a matrix-free action based on Fourier transforms. The numerical experiments reported here show that the preconditioner can reduce the condition number of the discrete (finite difference) Helmholtz operator by 10 orders of magnitude or more (see Tables 2-5) and improve the convergence of the GMRES iteration considerably (see Figures 1 and 3). The numerical experiments also show that the conditioning properties of the proposed method remain robust with respect to the points per wavelength and to the frequency ω𝜔\omegaitalic_ω. This latter observation is consistent with the pseudodifferential derivation whose neglected terms decay as the frequency increases. See remarks above Theorem 1. A second set of numerical experiments (see Subsection 5.2) illustrate the ability of the preconditioner to be implemented for wave scattering problems where an absorbing boundary layer/sponge is needed to truncate the physical domain and approximate the effect of the Sommerfeld radiation condition.

Finally we wish to list some limitations of the proposed preconditioning method and directions of possible improvement or extensions.

  1. 1.

    In general, the proposed pseudodifferential preconditioner is limited by its underlying technical assumptions, including the smoothness of the wave speed and other media properties. Although the principal symbol (21) does not contain derivatives of the acoustic properties, the next symbol does. In fact, we have observed a deterioration in the effectiveness of the preconditioner in Table 5 as the transition in wave speed values becomes sharper.

  2. 2.

    The domain of interest was discretized using a uniform Cartesian grid which is ideal for discretizing the Helmholtz operator using finite differences and for the implementation of the FFT for the application of the preconditioner. It would be interesting to investigate how to implement this preconditioner on non-Cartesian grids which are typically employed for finite element methods.

  3. 3.

    The speed of the proposed preconditioner relies on the application of the FFT. Unfortunately, FFT algorithms are not very well-suited for parallelization and distributed memory architectures. Hence, the FFT may prevent the proposed preconditioner from being fully parallelized.

  4. 4.

    Assumption (3), embedded in (25), may be too restrictive in some interesting cases. For instance, the wave speed c=c(x)𝑐𝑐𝑥c=c(x)italic_c = italic_c ( italic_x ) may be completely independent from the damping coefficient a=a(x)𝑎𝑎𝑥a=a(x)italic_a = italic_a ( italic_x ). In such cases, the interpolation method described in Section 4 should be extended to a bi-variate setting at the expense of increasing the number of interpolation points to sample the set [cmin,cmax]×[amin,amax]subscript𝑐minsubscript𝑐maxsubscript𝑎minsubscript𝑎max[c_{\rm min},c_{\rm max}]\times[a_{\rm min},a_{\rm max}][ italic_c start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_c start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ] × [ italic_a start_POSTSUBSCRIPT roman_min end_POSTSUBSCRIPT , italic_a start_POSTSUBSCRIPT roman_max end_POSTSUBSCRIPT ]. However, in many applications, such as ultrasonics in soft biological tissues, there is a one-to-one relationship between the wave speed and the attenuation coefficient of various tissues [43, 44]. Hence, the coefficient a𝑎aitalic_a can be expressed as a function of the wave speed c𝑐citalic_c.

  5. 5.

    In this paper, we have incorporated the presence of an absorbing layer/sponge defined by a complex-valued wave speed (see Subsections 4.4 and 5.2). However, perfectly matched layers (PML) have been shown to absorb the outgoing waves better. Hence, it would be interesting to extend the preconditioner to handle PMLs.

  6. 6.

    We have implemented the interpolation scheme using piecewise linear polynomials. Other interpolation bases such as piecewise Lagrange or Hermite polynomials could be implemented to improve the accuracy estimate of Theorem 2. In particular, Hermite interpolation may bring about interesting challenges as the derivative of a symbol with respect to the wave speed will introduce new symbols into the formulation.

  7. 7.

    It may be possible to design more accurate preconditioners by including more terms from the expansion (20)-(23). However, we anticipate some challenges. First, higher derivatives of the media properties (wave speed and damping coefficient) would need to be discretized. Second, a simple look at (22) reveals that the symbol q3subscript𝑞3q_{-3}italic_q start_POSTSUBSCRIPT - 3 end_POSTSUBSCRIPT represents a third order derivative in its numerator. This factor is eventually dominated by the denominator for very high frequencies |ξ|ω/cmuch-greater-than𝜉𝜔𝑐|\xi|\gg\omega/c| italic_ξ | ≫ italic_ω / italic_c. But for frequencies in the range 0|ξ|ω/cless-than-or-similar-to0𝜉less-than-or-similar-to𝜔𝑐0\lesssim|\xi|\lesssim\omega/c0 ≲ | italic_ξ | ≲ italic_ω / italic_c, we can expect the naive application of q3subscript𝑞3q_{-3}italic_q start_POSTSUBSCRIPT - 3 end_POSTSUBSCRIPT to behave unstably. Therefore, regularization or appropriate filters may be needed.

  8. 8.

    The Helmholtz equation is one of the simplest models for time-harmonic wave propagation. We believe that the work developed here could be extended to other models such as the Helmholtz equation is divergence form, Maxwell equations of electromagnetism, Navier equations of linear elasticity, and also to integral formulations such as the Lippmann-Schwinger equation.


Acknowledgment

The work of S. Acosta and T. Khajah was partially supported by NIH award 1R15EB035359-01A1. The work of B. Palacios was partially supported by Agencia Nacional de Investigación y Desarrollo (ANID) de Chile, Grant FONDECYT Iniciación N11220772. S. Acosta would like to thank the support and research-oriented environment provided by Texas Children’s Hospital.


References

  • [1] S. Acosta, J. Chan, R. Johnson, and B. Palacios. Pseudodifferential Models for Ultrasound Waves With Fractional Attenuation. SIAM J. Appl. Math., 84(4):1609–1630, 2024.
  • [2] H. Ammari, editor. Mathematical modeling in biomedical imaging I, volume 1983 of Lecture Notes in Mathematics. Springer-Verlag Berlin Heidelberg, 2009.
  • [3] H. Ammari, editor. Mathematical modeling in biomedical imaging II, volume 2035 of Lecture Notes in Mathematics. Springer-Verlag Berlin Heidelberg, 2012.
  • [4] D. Appelo, F. Garcia, and O. Runborg. WaveHoltz: Iterative solution of the Helmholtz equation via the wave equation. SIAM J. Sci. Comput., 42(4):A1950–A1983, 2020.
  • [5] K. Atkinson and W. Han. Theoretical Numerical Analysis: A Functional Analysis Framework, volume 39 of Texts in Applied Mathematics. New York : Springer-Verlag, 2001.
  • [6] J.-F. Aubry, O. Bates, C. Boehm, K. Butts Pauly, D. Christensen, C. Cueto, P. Gélat, L. Guasch, J. Jaros, Y. Jing, R. Jones, N. Li, P. Marty, H. Montanaro, E. Neufeld, S. Pichardo, G. Pinton, A. Pulkkinen, A. Stanziola, A. Thielscher, B. Treeby, and E. van’t Wout. Benchmark problems for transcranial ultrasound simulation: Intercomparison of compressional wave models. J. Acoust. Soc. Am., 152(2):1003–1019, 2022.
  • [7] G. Bao and W. Symes. Computation of pseudo-differential operators. SIAM J. Sci. Comput., 17(2):416–429, 1996.
  • [8] M. Bristeau, R. Glowinski, and J. Périaux. Controllability methods for the computation of time-periodic solutions; application to scattering. J. Comput. Phys., 147(2):265–292, 1998.
  • [9] Z. Chen and X. Xiang. A source transfer domain decomposition method for Helmholtz equations in unbounded domain. SIAM J. Numer. Anal., 51(4):2331–2356, 2013.
  • [10] P.-H. Cocquet and M. Gander. How large a shift is needed in the shifted Helmholtz preconditioner for its effective inversion by multigrid? SIAM J. Sci. Comput., 39(2):A438–A478, 2017.
  • [11] D. Colton and R. Kress. Inverse Acoustic and Electromagnetic Scattering Theory. Applied Mathematical Sciences. Springer, New York, NY, 3rd edition, 2013.
  • [12] L. Demanet and L. Ying. Discrete symbol calculus. SIAM Rev., 53(1):71–104, 2011.
  • [13] L. Demanet and L. Ying. Fast wave computation via Fourier integral operators. Math. Comput., 81(279):1455–1486, 2012.
  • [14] Pavel Drabek and Jaroslav Milota. Methods of Nonlinear Analysis. Birkhäuser Verlag AG, Basel, 2007.
  • [15] B. Engquist and L. Ying. Sweeping preconditioner for the Helmholtz equation: Hierarchical matrix representation. Commun. Pure Appl. Math., 64(5):697–735, 2011.
  • [16] B. Engquist and L. Ying. Sweeping preconditioner for the Helmholtz equation: Moving perfectly matched layers. Multiscale Model. Simul., 9(2):686–710, 2011.
  • [17] Y. Erlangga. Advances in iterative methods and preconditioners for the Helmholtz equation. Arch. Comput. Methods Eng., 15(1):37–66, 2008.
  • [18] Y. Erlangga, C. Vuik, and C. Oosterlee. On a class of preconditioners for solving the Helmholtz equation. Appl. Numer. Math., 50(3-4):409–425, 2004.
  • [19] O. Ernst and M. Gander. Why it is difficult to solve Helmholtz problems with classical iterative methods. Lect. Notes Comput. Sci. Eng., 83:325–363, 2012.
  • [20] M. Eslaminia and M. Guddati. A double-sweeping preconditioner for the Helmholtz equation. J. Comput. Phys., 314:800–823, 2016.
  • [21] M. Gander and F. Nataf. AILU: A preconditioner based on the analytic factorization of the elliptic operator. Numer. Linear Algebr. with Appl., 7(7-8):543–567, 2000.
  • [22] M. Gander and F. Nataf. An incomplete LU preconditioner for problems in acoustics. J. Comput. Acoust., 13(3):455–476, 2005.
  • [23] M. Gander and H. Zhang. A class of iterative solvers for the Helmholtz equation: Factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods. SIAM Rev., 61(1):3–76, 2019.
  • [24] M. Ganesh and C. Morgenstern. A coercive heterogeneous media Helmholtz model: formulation, wavenumber-explicit analysis, and preconditioned high-order FEM. Numer. Algorithms, 83(4):1441–1487, 2020.
  • [25] R. Glowinski, J.-L. Lions, and J. He. Exact and Approximate Controllability for Distributed Parameter Systems : A Numerical Approach, volume 117 of Encyclopedia of Mathematics and its Applications. Cambridge University Press, 2008.
  • [26] S. Gong, I. Graham, and E. Spence. Domain decomposition preconditioners for high-order discretizations of the heterogeneous Helmholtz equation. IMA J. Numer. Anal., 41(3):2139–2185, 2021.
  • [27] A. Grigis and J. Sjöstrand. Microlocal Analysis for Differential Operators: An Introduction. London Mathematical Society Lecture Note Series vol 196. Cambridge University Press, Cambridge, 1994.
  • [28] M. Grote and J. Tang. On controllability methods for the Helmholtz equation. J. Comput. Appl. Math., 358:306–326, 2019.
  • [29] E. Heikkola, S. Mönkölä, A. Pennanen, and T. Rossi. Controllability method for the Helmholtz equation with higher-order discretizations. J. Comput. Phys., 225(2):1553–1576, 2007.
  • [30] T. Kato. Perturbation Theory for Linear Operators. Springer Berlin Heidelberg, 2 edition, 1980.
  • [31] R. Kress. Numerical Analysis, volume 181 of Graduate Texts in Mathematics. Springer, New York, NY, 1998.
  • [32] R. Kress. Linear Integral Equations, volume 82 of Applied mathematical sciences. New York : Springer-Verlag, 2nd edition, 1999.
  • [33] W. McLean. Strongly Elliptic Systems and Boundary Integral Equations. Cambridge Univ. Press, 2000.
  • [34] A. Moiola and E. Spence. Is the Helmholtz equation really sign-indefinite? SIAM Rev., 56(2):274–312, 2014.
  • [35] J. Poulson, B. Engquist, S. Li, and L. Ying. A parallel sweeping preconditioner for heterogeneous 3D Helmholtz equations. SIAM J. Sci. Comput., 35(3):C194–C212, 2013.
  • [36] G. Soldati. Biomedical applications of ultrasound. In Compr. Biomed. Phys., volume 2, pages 401–436. Elsevier B.V., 2014.
  • [37] C. Stolk. A rapidly converging domain decomposition method for the Helmholtz equation. J. Comput. Phys., 241:240–252, 2013.
  • [38] C. Stolk. An improved sweeping domain decomposition preconditioner for the Helmholtz equation. Adv. Comput. Math., 43:45–76, 2017.
  • [39] C. Stolk. A time-domain preconditioner for the Helmholtz equation. SIAM J. Sci. Comput., 43(5):A3469–A3502, 2021.
  • [40] M. Taus, L. Zepeda-Núñez, R. Hewett, and L. Demanet. L-Sweeps: A scalable, parallel preconditioner for the high-frequency Helmholtz equation. J. Comput. Phys., 420:109706, 2020.
  • [41] M. E. Taylor. Pseudodifferential Operators and Nonlinear PDE. Progress in Mathematics vol 100. Birkhäuser, Boston, 1991.
  • [42] L. N. Trefethen and D. Bau. Numerical Linear Algebra. SIAM, 2022.
  • [43] T. Varslot and G. Taraldsen. Computer simulation of forward wave propagation in soft tissue. IEEE Trans. Ultrason. Ferroelectr. Freq. Control, 52(9):1473–1482, 2005.
  • [44] M. Verweij, B. Treeby, K. van Dongen, and L. Demi. Simulation of Ultrasound Fields. Elsevier B.V., 2014.
  • [45] L. Zepeda-Núñez and L. Demanet. The method of polarized traces for the 2D Helmholtz equation. J. Comput. Phys., 308:347–388, 2016.
  • [46] L. Zepeda-Núñez and L. Demanet. Nested domain decomposition with polarized traces for the 2D Helmholtz equation. SIAM J. Sci. Comput., 40(3):B942–B981, 2018.