Learning Dominant Wave Directions For Plane Wave Methods For High-Frequency Helmholtz Equations
Jun Fang, Jianliang Qian, Leonardo Zepeda-Núñez, Hongkai Zhao
LLearning Dominant Wave Directions For PlaneWave Methods For High-Frequency HelmholtzEquations
Jun Fang
Department of Mathematics, University of California, Irvine
Jianliang Qian
Department of Mathematics, Michigan State University
Leonardo Zepeda-N´u˜nez
Department of Mathematics, University of California, Irvine
Hongkai Zhao
Department of Mathematics, University of California, Irvine
Abstract
We present a ray-based finite element method (ray-FEM) by learning basis adap-tive to the underlying high-frequency Helmholtz equation in smooth media. Basedon the geometric optics ansatz of the wave field, we learn local dominant ray direc-tions by probing the medium using low-frequency waves with the same source. Oncelocal ray directions are extracted, they are incorporated into the finite element basisto solve the high-frequency Helmholtz equation. This process can be continued tofurther improve approximations for both local ray directions and the high frequencywave field iteratively. The method requires a fixed number of grid points per wave-length to represent the wave field and achieves an asymptotic convergence as thefrequency ω → ∞ without the pollution effect. A fast solver is developed for theresulting linear system with an empirical complexity O ( ω d ) up to a poly-logarithmicfactor. Numerical examples in 2D are presented to corroborate the claims. Keywords:
Helmholtz equation, numerical micro-local analysis, ray-FEM
Preprint submitted to Elsevier August 1, 2018 a r X i v : . [ m a t h . NA ] A ug . Introduction We consider the Helmholtz equation: H u := − ∆ u ( x ) − ω c ( x ) u ( x ) = f ( x ) , x ∈ Ω ⊆ R d , (1)plus boundary (or radiation) conditions, where ω is the frequency, c ( x ) > f ( x ) is the source distribution.The numerical solution of the Helmholtz equation (1) in the high-frequencyregime, i.e., ω (cid:29)
1, is notoriously hard to compute. From Shannon’s samplingprinciple [25], to resolve a general wave field oscillating at frequency ω , a meshsize h = O ( ω − ) is necessary and sufficient. Hence the intrinsic degrees of freedomis O ( ω d ), implying that the theoretical optimal overall complexity to solve (1) is O ( ω d ). In general, an overall complexity of optimal order is difficult to achieve inpractice due to two typical challenges: how to design a discretization method thatcan achieve both accuracy and stability with the minimum order of degrees of free-dom; and how to solve the linear system issued from the discretization with an almostlinear complexity as the frequency becomes large.Numerically, it is difficult to design a sparse discretization that can achieve bothaccuracy and stability under the condition ωh = O (1) as ω becomes large. This isusually called pollution effect in error estimates for finite element methods [2, 1, 14],i.e., the ratio between numerical error and best approximation error from a discretefinite element space is ω dependent.A popular strategy to solve the high-frequency Helmholtz equation (1) is to incor-porate appropriate oscillatory behaviors into the basis for Galerkin methods in orderto represent the oscillatory wave field governed by the Helmholtz equation more ac-curately and efficiently. The methods based on this strategy are usually called wavebased methods. The key issue for this strategy is how to design the oscillatory basis.Several approaches have been developed along this direction. One approach is toincorporate plane waves with a predetermined even distribution in directions intothe basis. For example, products of plane waves with local finite element basis areused in the generalized finite element method [1, 27]. Trefftz-type methods use localsolutions of the Helmholtz equation [12] as the basis functions, which in the case ofpiece-wise constant media are plane waves.A challenging issue in these approaches is how to choose the number of planewave directions `a priori . It is well known that in order to achieve a good accuracy,a fine, ω dependent, resolution in angle space is required. However, since manyof these directions may not be relevant to a specific problem setup, it will not only2ncrease the degrees of freedom significantly but also make the resulting linear systemill-conditioned due to the numerical loss of orthogonality of the elements of thebasis. For heterogeneous media, the construction of local solutions of the Helmholtzequation becomes more challenging, and some recent papers have focused on how toconstruct such local solution: a generalized plane wave method was developed in [15];and a two-scale method using local numerical solutions of the Helmholtz equationon a fine mesh as the basis function on a coarse mesh was developed in [9].Another approach to solve the high-frequency Helmholtz equation is based onthe geometric optics ansatz of the wave field. In the ansatz, the high-frequency wavefield is approximated by a superposition of a finite number, that we denote by N , ofdominant wave fronts of the form u ( x ) ≈ N (cid:88) n =1 A n ( x ) e iωφ n ( x ) (2)at each point; see [18] for examples. The main advantage is that both the amplitudefunctions A n ( x ) and phase functions φ n ( x ) are independent of the frequency ω andhence are non-oscillatory and smooth except at a measure zero set, e.g., focus points,caustics, corners in a smooth medium. However, both the amplitude and phasefunctions depend on the medium, source distribution and domain boundary, i.e.,they are global and problem specific. Once the phase functions of the wave fronts areavailable, the oscillatory pattern of the wave field is known and can be incorporatedinto the numerical approximation to improve both stability and accuracy. Instead ofincorporating many problem independent plane waves to the basis, the main idea ofphase based numerical methods is to explicitly incorporate the phase function intothe basis functions [10, 20].The main issue in phase-based methods is how to compute the phase functions,which are global unknown functions depending on the medium and the source dis-tribution. Only in simple geometry and homogenous medium the phase functioncan be computed analytically; otherwise, it has to be computed numerically [18].Usual approaches include Lagrangian methods, such as ray tracing by solving ordi-nary differential equations (ODE), and Eulerian methods based on partial differentialequation (PDE), e.g., the Eikonal equation, satisfied by the phase function. Boththese methods face numerical difficulties for general media in practice. For example,convergence and divergence of rays cause difficulties to obtain an accurate and uni-form ray field. Since the solution to the PDE is single valued, it can only capturethe phase function corresponding to the first arrival time of the wave front. For ageneral medium with varying speed, computing the appropriate global phase func-tions φ n ( x ) is a challenging task since different phase functions may be defined in3ifferent regions the boundaries of which are difficult to determine. In other words,the number of terms in the geometric optics ansatz varies from point to point; see[18] for examples. Moreover, since the phase functions are precomputed once for alland fed into the computation of the solution for the Helmholtz equation in theseapproaches, the phase function has to be computed extremely accurate because thenumerical error is amplified by ω . Instead of incorporating the global phase functioninto the basis, one can use local ray directions, (cid:98) d n ( x ) = ∇ φ n ( x ) |∇ φ n ( x ) | , as local dominantplane wave direction of e iω (cid:98) d n · x and incorporate them into the local basis. Comparedto the usual plane wave methods that typically require a large number of evenlydistributed wave directions independent of the problem, only directions of dominantwave fronts relevant to the problem are involved. Hence the degrees of freedom canbe kept minimal and ill-conditioning of the resulting linear system due to redun-dancy can be reduced. This combines the advantages of plane wave methods andphase-based methods. As is shown in [7], using dominant wave directions in planewave methods can significantly improve efficiency and accuracy for solving the high-frequency Helmholtz equation in heterogeneous medium. Moreover, it was shownthat there is some tolerance in the approximation of the dominant wave directions.For those examples in [7], a simple ray tracing method was used to find the rays andthe corresponding dominant wave directions at each point for a very simple setup ina homogeneous medium and the exact phase function was given for a heterogeneousmedium.Although incorporating problem specific dominant wave directions into basisfunctions has shown very promising results, a key issue in practice is how to findthe dominant wave directions in a stable, efficient and systematic way for generalmedia.In this work we propose an approach that is based on learning the dominant wavedirections specific to the medium and source distribution. In particular, we probe thesame medium by the same source, i.e., solving the Helmholtz equation (1) with thesame c ( x ) , f ( x ), for a relative low frequency (cid:101) ω ∼ √ ω . The computed wave field ispost-processed by numerical micro-local analysis (NMLA) or other signal processingtools to estimate dominant wave directions. These estimated wave directions are thenused in plane wave methods to solve the original high-frequency Helmholtz equation.In our approach, global phase functions are not needed. Instead, dominant wavedirections are computed locally where both the number of dominant wave directionsand the directions can vary from point to point. It provides the flexibility to dealwith general media. Moreover, once a more accurate wave field is computed, it canbe used to get a better estimation of the dominant wave directions and then usedto improve the high-frequency wave field again. In other words, both the dominant4ave directions and the high-frequency wave field can be computed and improved inan iterative way. In this paper, we develop a simple ray-based finite element (ray-FEM) method in 2D for smooth media as a proof of concept study of our proposedapproach. We start with a finite element mesh with mesh size h satisfying wh = O (1),i.e., a few points per wavelength. First, the low frequency Helmholtz equation with (cid:101) ω ∼ √ ω is solved on the mesh with quasi-optimality since (cid:101) ω h = O (1). Then NMLA[3, 5] (see Section 3) is applied to the computed low-frequency wave field to estimatelocal dominant wave directions. In the second step we incorporate these estimateddominant wave directions into the local finite element basis as the generalized finiteelement method to solve the high-frequency Helmholtz equation on the same mesh.If necessary, the solution for the high-frequency Helmholtz equation can also beprocessed by NMLA to improve the estimate of local dominant wave directions whichcan be used to further improve the high-frequency solution. We also develop aniterative fast solver for the discretized linear system based on polarized traces whichcan achieve O ( N ) complexity with a possible poly-logarithmic factor for smoothmedia, where N is the total number of unknowns. We use numerical examples to showthat our approach can achieve the following goals: (1) a stable discretization withdegrees of freedom of the minimal order O ( ω d ), (2) a fast solver with approximatelinear complexity for the discretized linear systems, (3) an asymptotic convergencerate of O ( ω − ) as ω → ∞ .Here is an outline of this paper. We first describe the ray-FEM using geometricoptics ansatz as the motivation and study its approximation property in Section2. In Section 3 we introduce the NMLA with its stability and local ray directionerror analyzed in Appendix A and B. Section 4 provides the full presentation of thenumerical algorithm whose complexity is given in Section 5. Numerical results arepresented in Section 6. Conclusion and future works are summarized in Section 7.
2. The Ray-FEM Method
In this section we describe the ray-FEM method for the Helmholtz equation andits rationale based on the geometric optics ansatz. We explain briefly the ansatz andapproximate it locally via a superposition of plane waves propagating in dominantdirections. These plane waves with their associated dominant directions are incorpo-rated into the finite element basis functions to improve both stability and accuracyfor the computation of the solution to the high-frequency Helmholtz equation.In this section we suppose that the dominant directions are known exactly. How-ever, in Section 3 we will describe how to learn the dominant wave directions byprobing the medium using low-frequency waves.5e use the following boundary value problem in 2D for the illustration of ourmethod, (cid:26) − ∆ u − k ( x ) u = f, in Ω , ∂u∂n + iβk ( x ) u = g, on ∂ Ω , (3)where Ω is a bounded Lipschitz domain in R , and k ( x ) = ω/c ( x ) is the inhomo-geneous wave number. (3) is usually refered to as the Helmholtz equation withimpedance boundary conditions. This equation was chosen in order to easily imposeanother types of boundary conditions by modifying the coefficient β . Specifically,the Dirichlet boundary condition corresponds to β = ∞ and the first order absorbingboundary condition to β = ±
1. Moreover, it is easy to extend (3) to incorporateabsorbing boundary conditions implemented via PML [6], as it will be performed inthe numerical experiments in Section 6.3.
The standard derivation of geometric optics ansatz is the use of WKJB ap-proximation [23, 16] (or the L¨uneberg-Kline expansion [17]) for the solution to theHelmholtz equation (1): u ( x ) ∼ e iωφ ( x ) ∞ (cid:88) (cid:96) =0 A (cid:96) ( x ) ω (cid:96) . (4)By taking ω → ∞ and considering only the first term one has u ( x ) = A ( x ) e iωφ ( x ) + O (cid:18) ω (cid:19) , (5)where A is usually called the amplitude and φ the phase. The key features of thegeometric optics ansatz are: • A and φ are independent of the frequency ω ; • A and φ depend on the medium, c ( x ); and the source distribution, f ( x ).Moreover, except for a small set of points, e.g., source/focus points, caustics, anddiscontinuities of the medium, they are smooth functions satisfying the followingPDE system:(eikonal) |∇ φ | = 1 c , (transport) 2 ∇ φ · ∇ A + A ∆ φ = 0 . (6)In general, the phase function, φ , and the amplitude function, A , are multi-valuedfunctions corresponding to multiple arrivals of wave fronts. Hence one can further6ecompose the geometric optics ansatz into a superposition of a number of wavefronts in the form: u ( x ) = N ( x ) (cid:88) n =1 A n ( x ) e iωφ n ( x ) + O (cid:18) ω (cid:19) , (7)where N ( x ) is the number of fronts/rays passing through x , and the phases φ n andamplitudes A n are single valued functions satisfying the eikonal/transport equations(6), each defined in a suitable domain with suitable boundary conditions [4].Based on the above geometric optics ansatz, one can derive a local plane waveapproximation at any point where φ n and A n are smooth with variations on a O (1)scale. Indeed, using Taylor expansions on a small neighborhood around an observa-tion point x , we have u ( x )= N ( x ) (cid:88) n =1 ( A n ( x )+ ∇ A n ( x )( x − x )) e iω ( φ n ( x )+ ∇ φ ( x ) · ( x − x )) + O (cid:18) h + ωh + 1 ω (cid:19) , (8)for | x − x | < h (cid:28) (cid:98) d n := ∇ φ n ( x ) |∇ φ n ( x ) | = c ( x ) ∇ φ n ( x ) (9)as the ray directions of the wave fronts at x , k ( x ) = ω/c ( x ), and B n ( x ) = ( A n ( x ) + ∇ A n ( x )( x − x )) e iω ( φ n ( x ) −∇ φ ( x ) · x ) (10)the affine complex amplitude. By replacing (9) and (10) in (8) we have u ( x ) = N ( x ) (cid:88) n =1 B n ( x ) e ik ( x ) (cid:98) d n · x + O (cid:18) h + ωh + 1 ω (cid:19) , (11)for | x − x | < h (cid:28)
1. From (11) we have that u can be approximated locally by asuperposition of plane waves propagating in certain directions with affine complexamplitudes. Moreover, as ω → ∞ , such that ωh = O (1), the asymptotic error forthe local plane wave approximation (11) is O ( ω − ), which is of the same order as theasymptotic error for the original geometric ansatz (7). We use (11) as the motivationto construct local finite element basis with mesh size h = O ( ω − ), in which an affinefunction is multiplied by plane waves oscillating in those ray directions, resulting inlocal approximations similar to (11). 7 .2. Ray-based FEM formulation We use a finite element method to compute the solution to (3) whose standardweak formulation is given byFind u ∈ H (Ω) , such that B ( u, v ) = F ( v ) , ∀ v ∈ H (Ω) , (12)where B ( u, v ) := (cid:90) Ω ∇ u · ∇ vdV − (cid:90) Ω k uvdV + iβ (cid:73) ∂ Ω kuvdS, (13) F ( v ) := (cid:90) Ω f vdV + (cid:73) ∂ Ω gvdS. (14)The domain, Ω, is discretized with a standard regular triangulated mesh, withmesh size h , which we denote by T h = { K } , where K represents a triangle of themesh. Using the aforementioned mesh we define two approximation spaces for thevariational formulation (12): • Standard FEM (S-FEM), we use low order P • Ray-FEM, we use P K ∈ T h , we denote by V j and x j , j = 1 , ,
3, the vertices of K and their coordinates respectively. Moreover, we denote by { ϕ j ( x ) } j =1 a partition ofunity consisting of piecewise bilinear functions satisfying ϕ j ( x i ) = δ ij , i, j = 1 , , δ ij is the Kronecker delta. The basis given by { ϕ j ( x ) } j =1 is usually called thenodal basis for Lagrange P V S ( K ) = span { ϕ j ( x ) , j = 1 , , } , (15)and the global P V S ( T h ) = { v ∈ C (Ω) : v | K ∈ V S ( K ) , ∀ K ∈ T h } . (16)To define the ray-FEM we enrich the P { (cid:98) d j,l } n j l =1 be n j ray directions at the vertex V j , define the ray-basedlocal approximation space by V Ray ( K ) = span { ϕ j ( x ) e ik j (cid:98) d j,l · x , k j = k ( x j ) , j = 1 , , , l = 1 , ..., n j } , and the global ray-FEM space by V Ray ( T h ) = { v ∈ C (Ω) : v | K ∈ V Ray ( K ) , ∀ K ∈ T h } .
8e can define the Standard FEM method byFind u ∈ V S ( T h ) , such that B ( u, v ) = F ( v ) , ∀ v ∈ V S ( T h ) . (17)Analogously, we define the ray-FEM method byFind u ∈ V Ray ( T h ) , such that B ( u, v ) = F ( v ) , ∀ v ∈ V Ray ( T h ) . (18) We provide a simple computation to estimate the approximation error of the ray-FEM space. In particular, we compute an asymptotic bound on inf u h ∈ V Ray ( T h ) || u − u h || L (Ω) , which is achieved by estimating the interpolating error using V Ray ( K ) as abasis.In the computation we assume that the ray direction, which is the gradient of thephase function φ , and the phase function itself, are exactly known. For simplicity,we assume N = 1 for the asymptotic formula in (7), i.e., that only one ray crosseseach point of the domain. In addition we suppose that f , the source, is zero insidethe domain. Under those circumstances A and φ are smooth; given that the source isoutside the domain there is no singularity in the amplitude, and given that only oneray crosses each point in the domain, no caustic occurs. From the geometric opticsansatz, we have u ( x ) = A ( x ) e iωφ ( x ) + O (cid:0) ω − (cid:1) . (19)Let Ω be our domain of interest, T h a triangle mesh of Ω with width h . Wedenote by N h the total number of vertices on the mesh T h , { x j } N h j =1 and { ϕ j ( x ) } N h j =1 are the coordinates of all mesh nodes and their corresponding nodal basis functionsfor standard P e iω [ φ ( x j ) −∇ φ ( x j ) · x j ] is a constant for the nodal basis associated to x j in an element K . From this observation we can easily deduce that the local ray-FEMspace can be rewritten as V Ray ( K ) = span { ϕ j ( x ) e ik j d j · x } = span { ϕ j ( x ) e iω ∇ φ ( x j ) · x } = span { ϕ j ( x ) e iω ∇ φ ( x j ) · x e iω [ φ ( x j ) −∇ φ ( x j ) · x j ] } = span { ϕ j ( x ) e iω [ φ ( x j )+ ∇ φ ( x j ) · ( x − x j )] } . Hence the nodal interpolation of the solution can be written as u I = N h (cid:88) j =1 A ( x j ) ϕ j ( x ) e iω [ φ ( x j )+ ∇ φ ( x j ) · ( x − x j )] , (20)9hich, by construction lies within the global ray-FEM space V Ray ( T h ).Let S j be the support of ϕ j ( x ), and | S j | ∼ O ( h ) be the area of S j . Then usingthe triangular inequality and the smoothness assumptions we have (cid:107) u − u I (cid:107) L (Ω) ≤ (cid:107) A ( x ) e iωφ ( x ) − (cid:80) N h j =1 A ( x j ) ϕ j ( x ) e iωφ ( x ) (cid:107) L (Ω) + (cid:107) (cid:80) N h j =1 A ( x j ) ϕ j ( x ) (cid:0) e iωφ ( x ) − e iω [ φ ( x j )+ ∇ φ ( x j ) · ( x − x j )] (cid:1) (cid:107) L (Ω) + O ( ω − ) ≤ (cid:107) A ( x ) − (cid:80) N h j =1 A ( x j ) ϕ j ( x ) (cid:107) L (Ω) + (cid:80) N h j =1 (cid:107) A (cid:107) L ∞ (Ω) (cid:107) e iωφ ( x ) − e iω [ φ ( x j )+ ∇ φ ( x j ) · ( x − x j )] (cid:107) L ( S j ) + O ( ω − ) (cid:46) h | A | H (Ω) + (cid:80) N h j =1 (cid:107) A (cid:107) L ∞ (Ω) ωh (cid:107)∇ φ (cid:107) L ∞ (Ω) | S j | + O ( ω − ) (cid:46) h | A | H (Ω) + ωh (cid:107) A (cid:107) L ∞ (Ω) (cid:107)∇ φ (cid:107) L ∞ (Ω) + O ( ω − ) . This impliesinf u h ∈ V Ray ( T h ) || u − u h || L (Ω) (cid:46) h | A | H (Ω) + ωh (cid:107) A (cid:107) L ∞ (Ω) (cid:107)∇ φ (cid:107) L ∞ (Ω) + O ( ω − ) . (21)Or, asymptotically, inf u h ∈ V Ray ( T h ) || u − u h || L (Ω) (cid:46) O ( h + ωh + ω − ) . (22)If the exact rays are known and the mesh size follows h ∼ ω − , then we haveinf u h ∈ V Ray ( T h ) || u − u h || L (Ω) (cid:46) O ( ω − ) , (23)i.e. that the approximation error decays linearly with ω , without oversampling. Remark 1.
The ray information can be incorporated into other Galerkin basis in thesame fashion. For example, in the hybrid numerical asymptotic method of [11], thebasis functions are constructed by multiplying nodal piece-wise bilinear function tooscillating functions with phase factors; the plane wave DG method of [7] employs theproducts of small degree polynomials and dominant plane waves as basis functions;the phase-based hybridizable DG method of [20] considers basis functions as productsof polynomials and phase-based oscillating functions. However, the phase or rayinformation in these methods is obtained from solving the eikonal equation with raytracing and related techniques. . Learning Local Dominant Ray Directions In Section 2 we use geometric optics to provide the motivation for the ray-FEM bybuilding an adaptive approximation space that incorporates ray information specificto the underlying Helmholtz equation. However, the ray directions, which dependon the medium and source distribution, are unknown quantities themselves, hencethey need to be computed or estimated. One way is to compute the global phasefunction, by either ray tracing or solving the eikonal equation, and take its gradient.As discussed in the introduction, computing the global phase function in a generalvarying medium can be difficult.In the present paper we propose a totally different approach. This novel approachis based on learning the dominant ray directions by probing the same medium withthe same source using a relative low-frequency wave. To be more specific, we firstsolve the Helmholtz equation (1) with the same speed function c ( x ), right handside f ( x ) and boundary conditions but with a relative low-frequency (cid:101) ω ∼ √ ω ona mesh with size h = O ( (cid:101) ω − ) = O ( ω − ) with a standard finite element method,which is quasi-optimal in that regime. Then the local dominant ray directions areestimated based on the computed low-frequency wave field. The key point is thatthe low-frequency wave has probed the medium specific to the problem globally andonly local dominant ray directions need to be learned, which allows us to handlemultiple arrivals of wave fronts seamlessly. In particular, we use numerical micro-local analysis (NMLA), which is simple and robust, in this work to extract thedominant ray directions locally. However, this is a signal processing task that can beaccomplished using other methods such as Prony’s method [8], Pisarenko’s method[22], MUSIC [24], matrix pencil [13], among many others. In this subsection, for the sake of completness, we provide a brief introduction toNMLA developed in [3, 5]. If we suppose that a wave field is a weighted superpositionof plane waves with the same wave number, but propagating in different directions,then the aim of NMLA is to extract from samples of the wave field, the directionsand the weights. In the sequel we use a 2D example to illustrate the method.Suppose that a wave field, denoted by u ( x ), is composed of N plane waves aroundan observation point x , u ( x ) = N (cid:88) n =1 B n e ik ( x − x ) · (cid:98) d n , | (cid:98) d n | = 1 . (24)11e suppose that we can sample the wave field, u ( x ), and its derivative on a circle S r ( x ) centered at x with radius r . The wave field can be written under the modelassumption in (24) as u ( x + r (cid:98) s ) = N (cid:88) n =1 B n e iα (cid:98) s · (cid:98) d n , α = kr, (cid:98) s ∈ S . (25)Furthermore, define implicitly the angle variable θ = θ ( (cid:98) s ); accordingly we denote θ n = θ ( (cid:98) d n ) the angles to recover, and x ( θ ) = x + r (cid:98) s ( θ ). Using the angle basednotation we sample the impedance quantity U ( θ ) = 1 ik ∂ r u ( x ( θ )) + u ( x ( θ )) , (26)which removes any possible ambiguity due to resonance, and improves the robustnessto noise for solutions to the Helmholtz equation [5], on the circle S r ( x ). Then weapply the filtering operator B to the impedance quantity B U ( θ ) := 12 L α + 1 L α (cid:88) l = − L α ( F U ) l e ilθ ( − i ) l ( J l ( α ) − iJ (cid:48) l ( α )) , (27)where L α = max(1 , [ α ] , [ α + ( α ) − . J l is the Bessel functions of order l , J (cid:48) l istheir derivatives and ( F U ) l = 12 π (cid:90) π U ( θ ) e − ilθ dθ (28)is the l -th Fourier coefficient of U . It is shown in [5] that B U ( θ ) = N (cid:88) n =1 B n S L α ( θ − θ n ) , (29)where S L ( θ ) = sin([2 L +1] θ/ L +1] sin( θ/ . As a consequence, we have that if α = kr → ∞ thenlim α →∞ B U ( θ ) = (cid:26) B n , if θ = θ n (or (cid:98) s = (cid:98) d n );0 , otherwise . (30)Then it is possible to obtain the directions and the amplitudes by picking the peaksin the filtered data in (29).However, for applications, the measured data is never a perfect superposition ofplane waves; therefore, we provide, for completeness, stability and error estimates for12MLA from [5] in Appendix A. In principle, as long as the perturbation is relativesmall with respect to the true signal, the estimation error is O ( kr ). In other words,the larger the radius of the circle compared to wavelength the more accurate theestimation is.In our application, our data is the numerical solution of the Helmholtz equation.In addition to noises and numerical errors, there are perturbations due to two modelerrors : • the geometric optics ansatz has an asymptotic error of order O ( ω − ) (see (7)); • in the geometric ansatz the wave field at a point is a superposition of curvedwave fronts. In particular, the curvature of the wave fronts results in a com-promise in the choice of the radius of the sampling circle to be of order O ( ω − )for the NMLA in order to achieve the stability and the minimal error of order O ( ω − ).Detailed analysis is provided in Appendix B. Below is a summary of the NMLAalgorithm. Algorithm 1
NMLA function d ω = NMLA ( x , T h , ω, h, c, u ) choose r ∼ ω − (cid:46) Radius for the sampling circle choose M ∼ ωr (cid:46) Number of sampling points ∆ θ = 2 π/M , (cid:46) Angular discretization for θ = 0 : ∆ θ : 2 π do x ( θ ) = x + r (cid:98) s ( θ ) U ( θ ) = ωic ( x ) ∂ r u ( x ( θ )) + u ( x ( θ )) (cid:46) Sample impedance data end for F ( θ ) = B U ( θ ) (cid:46) Apply the filter (27) θ est = SharpPeakLocations ( θ, F ( θ )) d ω = d ( θ est ) end function In this section we incorporate the errors in the estimation of the ray directionsinto the approximation error for the ray-FEM method, in which ray directions arefirst estimated by using Algorithm 1 to the solution of the Helmholtz equation witha relative low-frequency. The estimated ray directions are then used to generate the13pproximation space. With the same assumptions as in Section 2.3, we estimate anupper bound on inf u h ∈ V hRay ( T h ) (cid:107) u − u h (cid:107) L (Ω) , (31)when the ray-FEM space, V hRay ( T h ), is constructed using the estimated ray directionsfrom high-frequency waves by NMLA.From Appendix B, the error estimation of dominant ray directions is O ( ω − / ).The numerical ray-FEM space V hRay ( T h ) is defined similarly as V Ray ( T h ) with the exactray directions { (cid:98) d j } replaced by the ones { (cid:98) d hj } estimated by NMLA and | (cid:98) d j − (cid:98) d hj | ∼O ( ω − / ).We denote by u hI = N h (cid:88) j =1 A ( x j ) ϕ j ( x ) e iω [ φ ( x j )+1 /c ( x j ) (cid:98) d hj · ( x − x j )] (32)the nodal interpolation of the solution in V hRay ( T h ). Then we have (cid:107) u I − u hI (cid:107) L (Ω) = (cid:107) (cid:80) N h j =1 A ( x j ) ϕ j ( x ) e iωφ ( x j ) ( e iω ∇ φ ( x j ) · ( x − x j ) − e iω/c ( x j ) (cid:98) d hj · ( x − x j ) ) (cid:107) L (Ω) ≤ (cid:80) N h j =1 (cid:107) A (cid:107) L ∞ (Ω) (cid:107) e iω/c ( x j ) (cid:98) d j · ( x − x j ) − e iω/c ( x j ) (cid:98) d hj · ( x − x j ) (cid:107) L ( S j ) (cid:46) (cid:80) N h j =1 (cid:107) A (cid:107) L ∞ (Ω) ωh (cid:107) c − (cid:107) L ∞ (Ω) | (cid:98) d j − (cid:98) d hj || S j | (cid:46) ω / h (cid:107) A (cid:107) L ∞ (Ω) (cid:107) c − (cid:107) L ∞ (Ω) . Hence,inf u h ∈ V hRay ( T h ) (cid:107) u − u h (cid:107) L (Ω) ≤ (cid:107) u − u hI (cid:107) L (Ω) ≤ (cid:107) u − u I (cid:107) L (Ω) + (cid:107) u I − u hI (cid:107) L (Ω) (cid:46) h | A | H (Ω) + ωh (cid:107) A (cid:107) L ∞ (Ω) (cid:107)∇ φ (cid:107) L ∞ (Ω) + ω / h (cid:107) A (cid:107) L ∞ (Ω) (cid:107) c − (cid:107) L ∞ (Ω) + O ( ω − ) . (33)Or, more compactly, we have thatinf u h ∈ V hRay ( T h ) (cid:107) u − u h (cid:107) L (Ω) = O ( h + ωh + ω / h + ω − ) . (34)Comparing with (22) and (34), the error in the estimation of dominant ray di-rections due to NMLA leads to the extra term ω / h , which is the leading order inthe high-frequency regime. Specifically, if ωh = O (1), then we haveinf u h ∈ V hRay ( T h ) (cid:107) u − u h (cid:107) L (Ω) = O ( ω − / ) . (35)14 . Algorithms In this section we provide the full algorithm for ray-FEM including a fast iterativesolver based on a modification of the method of polarized traces for the resultinglinear systems. In order to streamline the presentation and to make the algorithmeasier to understand, we introduce several subroutines, which the main algorithmrelies on.We can separate the full algorithm into three conceptual stages:1. probing the medium by solving a relative low-frequency Helmholtz equationwith the standard FEM,2. learning the dominant ray directions from low-frequency probing wave field byNMLA,3. solving the high-frequency Helmholtz equation in ray-FEM space.If necessary the second stage can be applied to the high-frequency wave fieldcomputed in stage 3 to improve the estimation of dominant ray directions and thenrepeat stage 3.We remind the reader that the ultimate objective of the algorithm presented inthis paper (i.e., Algorithm 5) is to solve the Helmholtz equation (1) at frequency ω with a total O ( ω d ) (up to poly-logarithmic factors) computational complexity. Inorder to achieve this objective, we discretize the PDE with a mesh size h = O ( ω − ),which leads to a total of O ( ω d ) number of degrees of freedom and a sparse linearsystem with O ( ω d ) number of nonzeros. Then we develop a fast iterative solver withquasi-linear complexity to solve the resulting linear system after discretization. Belowis a more detailed description of the three stages. Finally, following the notationdefined in the prequel, we denote the triangular mesh by T h . We first solve the Helmholtz equation with frequency (cid:101) ω ∼ √ ω in the same mediumand with the same source on T h . The low-frequency problem is solved using the stan-dard finite element method (S-FEM) with linear elements as prescribed by Algorithm2. 15 lgorithm 2 Standard FEM Helmholtz Solver function u ω,h = S-FEM ( ω, h, c, f, g ) for i, j = 1 : N h doH i,j = B ( ϕ i , ϕ j ) (cid:46) Assemble Helmholtz matrix b j = F ( ϕ j ) (cid:46) Assemble right-hand side end foru ω,h = H − b (cid:46) Solve linear system end function
Let u (cid:101) ω,h = S-FEM ( (cid:101) ω, h, c, f, g ) denote the S-FEM solution of the low-frequencyHelmholtz equation on T h . Since ˜ ω h = O (1), S-FEM is quasi-optimal in the norm (cid:107) · (cid:107) H := (cid:107)∇ · (cid:107) L + k (cid:107) · (cid:107) L [19], and it has optimal L error estimate [28]. Once the low-frequency problem has been solved, we extract the dominant raydirections from u (cid:101) ω,h using NMLA as described in Section 3.1 around each mesh node.We utilize the smoothness of the phase functions, and hence the smoothness of theray directions field to reduce the computational cost. The reduction is achieved byrestricting the learning of the dominant ray directions to vertices of a coarse meshdown-sampled from T h . Such coarse mesh is denoted by T h c , where h c = O ( √ h ).The resulting dominant ray directions are then interpolated onto the fine mesh T h .Note that at each vertex of T h c , the wave field u (cid:101) ω,h on the fine mesh T h is usedfor NMLA to estimate the dominant ray directions. There are three sources of errorin the learning stage we need to be aware of: • numerical error of u (cid:101) ω,h , • model error in geometric optics ansatz, • interpolation error.The numerical error for u (cid:101) ω,h by Algorithm 2 in L norm [28] is O ( (cid:101) ωh + (cid:101) ω h ) = O ( ω − ), which is negligible with respect to the model error in the geometric opticsansatz. The error introduced by the geometric optics approximation is O ( (cid:101) ω − ) asshown in Section 3.1 and Appendix B. The error due to linear interpolation on T h c toget ray direction estimation at every vertex on T h is O ( h c ) = O ( h ) = O ( ω − ), whichis much smaller than the model error in geometric optics ansatz. Hence the overallerror in ray direction estimation based on NMLA on u (cid:101) ω,h and interpolation is O ( (cid:101) ω − ).The dominant ray direction estimation algorithm is summarized in Algorithm 3. For16ach node x j on mesh T h the number of dominant ray directions is denoted by n j , d jω,h = { d j,lω,h } n j l =1 . Algorithm 3
Ray Learning function { d jω,h } N h j =1 = RayLearning ( ω, h, h c , c, u ω,h ) for j = 1 : N h c do d jω,h c = NMLA ( x cj , T h , ω, h, c, u ω,h ) end for { d jω,h } N h j =1 = Interpolation ( T h c , T h , { d jω,h c } N hc j =1 ) end function Once the dominant ray directions on T h have been computed, we can construct theray-FEM space V Ray ( T h ) and solve the high-frequency Helmholtz equations following(18), which is implemented in Algorithm 4. Algorithm 4
Ray-FEM Helmholtz Solver function u d ω ,h = Ray-FEM ( ω, h, c, f, g, { d jω,h } N h j =1 ) N dof = 0 for j = 1 : N h , l = 1 : n j do N dof = N dof + 1, m = N dof ψ m ( x ) = ϕ j ( x ) e iw/c ( x j ) d j,l · x (cid:46) Construct ray-FEM basis functions (cid:98) ψ = ψ m ( x j ) (cid:46) Nodal values of ray-FEM basis functions end for for m, n = 1 : N dof do H m,n = B ( ψ m , ψ n ) (cid:46) Assemble Helmholtz matrix b n = F ( ψ n ) (cid:46) Assemble right-hand side end for v = H − b (cid:46) Coefficients of ray-FEM basis functions u d ω ,h = v · (cid:98) ψ (cid:46) Ray-FEM solution on mesh nodes end function
In general, the more accurate of the dominant ray directions learned by NMLA,the more accurate ray-FEM high-frequency solution can be computed by Algorithm4. From section 4.2, the accuracy order of learning stage from low-frequency wavefield is O ( (cid:101) ω − ), and following the error analysis of section 3.2, the consequent ray-FEM solution has the same order of accuracy. However, we can apply the learning17tage to the computed high-frequency wave field to improve approximation for domi-nant ray directions. Moreover, the improved ray information can be used to improvethe approximation for high-frequency wave field computed by ray-FEM. This processcan be continued iteratively to further improve the approximations of the ray direc-tions and the high-frequency wave field. Therefore, we develop an iterative ray-FEMHelmholtz solver for high-frequency ω in Algorithm 5. Algorithm 5
Iterative Ray-FEM High-Frequency Helmholtz Solver function u d ω ,h = IterRay-FEM ( ω, c, f, g ) (cid:101) ω ∼ √ ω , h ∼ ω − , h c ∼ ω − u (cid:101) ω,h = S-FEM ( (cid:101) ω, h, c, f, g ) (cid:46) Low-frequency waves { d (cid:101) ω,h } = RayLearning ( (cid:101) ω, h, h c , c, u (cid:101) ω,h ) (cid:46) Low-frequency ray learning u d (cid:101) ω ,h = Ray-FEM ( ω, h, c, f, g, { d (cid:101) ω,h } ) (cid:46) High-frequency waves tol = 1, niter = 0, u ω,h = u d (cid:101) ω ,h while tol > (cid:15) or niter > max iter do { d ω,h } = RayLearning ( ω, h, h c , c, u ω,h ) (cid:46) High-frequency ray learning u d ω ,h = Ray-FEM ( ω, h, c, f, g, { d ω,h } ), u ω,h = u d ω ,h tol = (cid:107) u ω,h − u ω,h (cid:107) L (Ω) / (cid:107) u ω,h (cid:107) L (Ω) niter = niter + 1, u ω,h = u ω,h end while end functionRemark 2. Extensive numerical experiments and Appendix A suggest that theNMLA process in learning dominant ray directions is remarkably stable even fornoisy plane wave data. Hence, the iterative process in Algorithm 5 usually needsvery few iterations to reach the desired accuracy. Typically, we only need or iterations in our numerical tests.4.4. Fast linear solver To achieve the overall complexity mentioned in the introduction, it is necessaryto solve the linear system resulting from both standard finite element methods andray-FEM, which we write in a generic form as Hu = f , (36)in linear complexity (up to poly-logarithmic factors). This solver is, in fact, thebottleneck of Algorithms 2 and 4. 18or a smooth medium, this can be achieved by modifying the method of polarizedtraces [30], of which we provide a brief review here. For further details we refer theinterested readers to [30]. The method of polarized traces is a domain decompositionmethod that encompasses the following aspects • layered domain decomposition; • absorbing boundary conditions between subdomains implemented via PML [6]; • transmission conditions issued from a discrete Green’s representation formula; • efficient preconditioner arising from localization of the waves via an incompleteGreen’s formula.The first two aspects can be effortlessly implemented. Consider a layered partitionof Ω into L slabs, or layers { Ω (cid:96) } L(cid:96) =1 . Define f (cid:96) as the restriction of f to Ω (cid:96) , i.e., f (cid:96) = f χ Ω (cid:96) ; and define the local Helmholtz operators as H (cid:96) u := (cid:0) −(cid:52) − mω (cid:1) u in Ω (cid:96) , (37)with absorbing boundary conditions implemented via PML between slabs.The method of polarized traces aims to solve the global linear system in (36) bysolving the local systems H (cid:96) , which are the discrete version of (37).In order to solve the global system, or in this case to find a good approximatesolution, we need to “glue” the subdomains together, which is achieved via a discreteGreen’s integral formula deduced by imposing discontinuous solutions.In the original formulation of the method of polarized traces [30], the Green’srepresentation formula was used to build a global surface integral equation (SIE) atthe interfaces between slabs. The SIE was solved using an efficient preconditionercoupled with a multilevel compression of the discrete kernels to accelerate the on-linestage of the algorithm. The original algorithm had superlinear off-line complexitywhich was amortized among a large number of right-hand sides, which represents atypical situation in explorations geophysics.In the context of the present paper, the linear systems issued from the ray-based FEM depends on the source distribution, making it impossible to amortizea super-linear off-line cost. In order to reduce the off-line cost we use a matrix-free formulation (see Chapter 2 in [29]) with a domain decomposition in thin layers.In this case, the cost per iteration is linear on the number of degrees of freedom,depending on the grow or the auxiliary degrees of freedom corresponding to thePML’s. Finally the convergence is normally achieved in O (log ω ) iterations, as itwill be shown in the sequel. 19 . Complexity In this section we provide an overall computational complexity count of our al-gorithm for the high-frequency Helmholtz equation (1) in terms of ω . Our countincludes learning ray directions by NMLA and the linear solver for the discretizedsystems from both standard FEM and ray FEM for low-frequency and high-frequencyHelmholtz equation respectively. The summary of the complexity of the method isgiven in Table 2. As described in Section 4.2, Algorithm 3 applies NMLA to computed wave fieldswith low-frequency (cid:101) ω ∼ √ ω or high-frequency ω and first estimate ray directions atvertices on a downsampled coarse mesh T h c and then interpolate the ray directionsto the vertices on a fine mesh T h . We remind the reader the following scalings: h = O ( ω − ), h c = O ( √ h ) = O ( ω − ) . These scalings allow us to strike a balanceamong the number of observation points at which NMLA is used to estimate raydirections, the radius of sampling circle, and the corresponding number of samplingpoints on the circle to resolve the wave field to reach the optimal accuracy of NMLAwith desired total computational complexity.We estimate the dominant computational cost which is applying NMLA to thehigh-frequency wave field to get ray direction estimation on T h in 2D as an illustrationof Table 1. It is shown in Appendix B, the least error that can be achieved byNMLA is O ( ω − ) when radius r of the sampling circle centered at an observationpoint is O ( ω − ). Hence the number of points sampled on the circle to resolve thewave field of frequency ω is M ω = O ( ω ). Since NMLA is a linear filter based ona Fourier transform in angle space, the corresponding computational complexity is O ( M ω log M ω ) [3]. The number of observation points we need to perform NMLAis the number of vertices on the coarse mesh which is O ( h − c ) = O ( ω ). Hence thecomputation cost to obtain ray directions at all vertices on the coarse mesh by NMLAis O ( ω log ω ). Finally the ray directions estimated at the vertices on the coarse meshby NMLA are interpolated to the fine mesh T h . Interpolation is a linear operationand hence its computation complexity is O ( ω ).Table 1 provides the complexity for each operation of NMLA, where d is thedimension and C NMLA , C ray,h c , C Int , and C ray,h is the computation complexity ofNMLA at a single vertex, NMLA on the under sampled coarse mesh, interpolationof local ray directions to the fine mesh and full algorithm for learning local raydirections at frequency ω on mesh T h respectively.20requency r M ω C NMLA C ray,h c C Int C ray,h ω ω − ω d − ω d − log ω ω d − log ω ω d ω d Table 1: Computational complexities of estimating ray directions on a mesh T h with h = O ( ω − ). The most computationally intensive component in the whole ray-FEM algorithmis solving the linear systems after discretization of the Helmholtz equation. Algorithm5 solves both u (cid:101) ω,h = S-FEM ( (cid:101) ω, h, c, f, g ) by S-FEM and u d ω ,h = Ray-FEM ( ω, h, c, f, g, { d jω,h } N h j =1 )by ray-FEM on the same mesh T h . Each solver is composed of three steps: the as-sembling step, the setup step, and the iterative solve step.Since the basis functions are locally supported, the resulting matrix is sparse.The complexity of the assembling step is of the same order as the degrees of freedom N h = O ( ω d ).In the setup stage, the computational domain is decomposed into subdomains ofthin layers whose width is comparable to the characteristic wavelength. The localproblems in each subdomain are factorized using multifrontal method in O ( √ N h )time (or O ( √ N h log N h ) time depending on the width of the auxiliary PML for eachsubdomain in term of the wavelength). Given that the layers are O (1) elementsthick, we have to factorize O ( √ N h ) subsystems, which results in a total O ( N h ) (or O ( N h log N h )) asymptotic complexity for the setup step.Finally, for the iterative solve step, each application of the preconditioner involves6 local solves per layer, each one performed with O ( √ N h ) ( or O ( √ N h log N h )) com-plexity. Then given that we have O ( √ N h ) layers, then we have an overall O ( N h )(or O ( N h log N h )) complexity per iteration. Extensive numerical experiments sug-gest that the number of iterations to converge is O (log N h ) for the both high- andlow-frequency solves for smooth media. Hence, the empirical overall complexity is O ( N h log N h ) for the high-frequency solve and O ( N h log N h ) for the low-frequencyone, which is summarized in the following table 2:Methods S-FEM Learning Ray-FEM Iterative Ray-FEMFrequency √ ω √ ω or ω ω ω Complexity O ( ω d log ω ) O ( ω d ) O ( ω d log ω ) O ( ω d log ω ) Table 2: Overall computational complexities with mesh size h = O ( ω − ). . Numerical Experiments In this section we provide several numerical experiments to test the proposedray-FEM and corroborate our claims. For all cases, the domain of interest is Ω =[ − / , / with different source terms and boundary conditions. Ω is discretizedusing a standard traingular mesh. The mass and stiffness matrices are assembledusing a high-order Gaussian quadrature rule to compute the integrals numerically . In the first test, the exact solution to the Helmholtz equation with Robin bound-ary condition is the wave field (normalized by the frequency ω ) corresponding to apoint source outside the domain. It is given by u ex ( x, y ) = √ ωH (1)0 ( ω (cid:112) ( x − + ( y − ) . (38)Numerically we solve the Helmholtz equation (1) with c ( x ) ≡
1, source f ( x ) ≡ ω and test convergence for both the raydirection estimation by NMLA and the final numerical solution by ray-FEM.First, a probing wave with low-frequency (cid:101) ω = √ ω is solved by standard FEM.Then NMLA is applied to the low-frequency probing wave to get an estimation of thelocal dominant ray directions d (cid:101) w . Instead of using the regular NMLA for plane wavedecomposition, we use NMLA with a curvature correction version [5] to estimate thenormal direction of a circular wave front. The local ray direction information is usedin ray-FEM to produce the first numerical solution to the high frequency Helmholtzequation u d (cid:101) w .We employ one more iteration in the framework of iterative ray-FEM by applyingNMLA to u d (cid:101) w to get an improved local ray direction estimation d w and then use itagain in ray-FEM to get a more accurate numerical solution to the high-frequencyHelmholtz equation u d w .Due to the curvature correction in NMLA, it can be shown that the error inray direction estimation (cid:107) d (cid:101) w − d ex (cid:107) ( (cid:107) d w − d ex (cid:107) ) is O ( (cid:101) ω − ) ( O ( ω − )) by using ananalysis similar to that in [5] and Appendix B. Using the estimate in Section 3.2,one can show that the approximation error for numerical ray-FEM space is at least O ( (cid:101) ω − ) ( O ( ω − )) if d (cid:101) w ( d w ) is used. Given the expression of the mass and stiffness matrices, which are polynomials times a planewave, it is possible to compute the integral analytically [21]. ωh = O (1). Quasi-optimality holds if the ratio ofnumerical solution error to the best approximation error is bounded by a constantwhich is independent of frequency ω . Here we assume the best approximation error inray-FEM space inf u h ∈ V hRay ( T h ) (cid:107) u − u h (cid:107) is of the same order as the interpolation error (cid:107) u − u hI (cid:107) . So we can define the estimated quasi-optimality constant by (cid:107) u − u h (cid:107) L (cid:107) u − u hI (cid:107) L ,and it is shown in Figure 2. Also it shows that one more iteration using iterative ray-FEM can significantly improve the final numerical solution to the order of O ( ω − ),which is of the same order when exact ray direction d ex is used in ray-FEM, due to theasymptotic error for geometric ansatz. Note, however, that NMLA with curvaturecorrection is only valid for a single point source in homogeneous media. ω/ π
20 40 80 1601 /h
120 240 480 960 (cid:107) θ ( d (cid:101) ω ) − θ ex (cid:107) L (cid:107) θ ( d ω ) − θ ex (cid:107) L (cid:107) u d (cid:101) ω − u ex (cid:107) L (cid:107) u d ω − u ex (cid:107) L (cid:107) u d ex − u ex (cid:107) L Table 3: Errors of one point source problem for fixed NPW = 6. θ ex is the exact ray angle, θ ( d (cid:101) ω )and θ ( d ω ) are ray angle estimations using low and high frequency waves, respectively; u d (cid:101) ω , u d ω and u d ex are ray-FEM solutions using low-frequency ray estimation d (cid:101) ω , high frequency ray estimation d ω , and exact ray d ex , respectively. Next we show that our method can handle multiple wave fronts by probing thewhole domain and extracting dominant ray directions locally. The setup is exactlyas above except that there are four point sources. The exact solution is given by u ex ( x, y ) = √ ωH (1)0 ( ω (cid:112) ( x + 20) + ( y + 20) ) + 2 √ ωH (1)0 ( ω (cid:112) ( x − + ( y − )+0 . √ ωH (1)0 ( ω (cid:112) ( x + 20) + ( y − ) − √ ωH (1)0 ( ω (cid:112) ( x − + ( y + 20) ) . (39)The main difficulty of this example compared to the one above, is that the low-frequency wave solution by the standard FEM contains multiple wave fronts at each23oint due to the interference of multiple sources. The numerical results are shownin right column of Figure 1. In this case, NMLA with curvature correction doesnot apply so we use the the standard NMLA version for plane wave decompositiondescribed in Section 3.1 to estimate local dominant ray directions. As analyzedin Section 3.2 and Appendix B, the expected error for ray direction estimationand numerical solution is of order O ( ω − / ) due to the curved wave fronts. Thenumerical results show that the ray-FEM meets the expectation without pollutionas the frequency increases. Here we show that by incorporating the estimated ray directions, ray-FEM cancapture the phase much more accurately. We test our algorithm with a source insidethe domain. In particular, we use a point source, given its importance in manypractical applications, in particular, in geophysics, in which the sources are oftenmodeled as point sources. Moreover, in applications oriented towards inverse andimaging problems, having a numerical method that produces the correct phase inthe far field is of great importance in order to properly locate features in the image.In this experiments we focus our attention on the far-field since our currentmethod can not deal with singularities in amplitude and phase at source points.We start by solving the Helmholtz equation using a slight modification of Algorithm5 and using standard finite elements at the source point. In this case the source termis located inside the domain f = δ ( x − ( − . , − . h ) as the right hand side. Forvertices near the source, we use the exact ray direction, the radial direction, in ourray-FEM. And we find the ray directions by NMLA for vertices away from the source;see Figure 3 left part for the ray direction field.To demonstrate the phase errors in numerical solutions, we plot the computedwave field on a 90 degree part of an annulus [26], with the radial coordinate varyingon an interval of about two wavelengths; see Figure 3 right part.In this case the frequency ω = 250 π is fixed, but we increase the number of gridpoints per wavelength. Figure 4 depicts the behavior of both ray-FEM solution andstandard FEM solution. From the figure we can easily observe the superiority of theray-FEM on minimizing the phase error, even using relatively coarse meshes.We then solve the Helmholtz equation using a heterogeneous medium given byFigure 5 with the source located inside. We also provide an experiment where weshow the ability of the method in this paper to handle wave field with caustics; seeFigure 6. Again radial directions are used for local ray directions near the sourcepoint. 24 .3. Complexity tests In this subsection we test the computational complexity for ray-FEM. A key stepis solving the sparse linear systems generated by ray-FEM using iterative methodswith a performant preconditioner, e.g., domain decomposition techniques coupledwith high-quality absorbing/transmission boundary conditions. In our tests, we usea modification of the method of polarized traces to solve the linear systems resultingfrom both standard FEM and ray-FEM as described in Section 4.4.We use the numerical experiments to demonstrate the overall computational com-plexity of our ray-FEM method. In particular, we solve the Helmholtz equation witha point source in both homogeneous medium and heterogeneous medium. We com-pute for many different frequencies, using Algorithm 5 with only one iteration of ray-FEM, the solution to the Helmholtz equation posed on Ω = [ − . , . × [ − . , . low = [ − , × [ − , O ( N ) up topoly-logarithmic factors as shown in our complexity study. The low-frequency solverhas a slightly higher asymptotic cost in this case, given the ratio between the widthof the PML and the characteristic wavelength inside the domain.Figure 8 shows the runtime for solving the Helmholtz equation with a point sourceinside a heterogenous medium. We can observe the same scaling as before, albeitwith slightly larger constants. 25 ω -4 E rr o r k θ ( d f ω ) − θ ex k L (Ω) k θ ( d ω ) − θ ex k L (Ω) O ( ω − ) ω -2 E rr o r k θ ( d f ω ) − θ ex k L (Ω) k θ ( d ω ) − θ ex k L (Ω) O ( ω − / ) ω -5 E rr o r k u d f ω − u ex k L (Ω) k u d ω − u ex k L (Ω) O ( ω − ) ω -4 E rr o r k u d f ω − u ex k L (Ω) k u d ω − u ex k L (Ω) O ( ω − / ) ω -5 E rr o r k u d ex − u ex k L (Ω) O ( ω − ) ω -6 E rr o r k u d ex − u ex k L (Ω) O ( ω − ) Figure 1: Tests with source outside domain, NPW = 6. Left: one point source; Right: fourpoint sources. Top: ray direction errors; Middle: errors of ray-FEM solutions with ray directionsestimated by NMLA; Bottom: errors of ray-FEM solutions with exact ray directions.
20 40 60 80 100 120 140 160 180 200
Frequency ω /2 π O p t i m a li t y r e l a t i on C op t Figure 2: The stars defined by (cid:107) u − u h (cid:107) L (cid:107) u − u hI (cid:107) L with NPW= 6, give an indication of the optimalityconstant. -0.6 -0.4 -0.2 0 0.2 0.4-0.6-0.4-0.200.20.4 Figure 3: One point source inside homogeneous medium domain, ω = 80 π , NPW = 6. Left: raydirection field captured by NMLA; Right: polar plot of ray-FEM solution, r/λ : the number ofwavelength away from the source. igure 4: Polar plot of ray-FEM solution ω = 250 π . r/λ : the number of wavelength away from thesource. Figure 5: One point source inside heterogeneous medium domain with Gaussian wave speed c ( x, y ) = 3 − . e − (( x +0 . +( y − . ) / . , ω = 80 π , NPW = 10. Left: ray direction field cap-tured by NMLA; Right: wave field computed by ray-FEM.Figure 6: One point source inside heterogeneous medium domain with sinusoidal wave speed c ( x, y ) = 1 + 0 . πx ), ω = 80 π , NPW = 10. Left: wave speed; Right: wave field computed byray-FEM. N = n T i m e [ s ] SolveSetup O ( N log N ) O ( N log N ) N = n T i m e [ s ] SolveSetup O ( N log N ) O ( N ) Figure 7: Runtime for solving the Helmholtz equation with homogeneous wave-speed using GMRESpreconditioned with the the method of polarized traces. The tolerance was set up to 10 − . Left:runtime for solving the low-frequency problem. Right: Runtime for solving the high-frequencyproblem with the adaptive basis. N = n T i m e [ s ] SolveSetup O ( N log N ) O ( N log N ) N = n T i m e [ s ] SolveSetup O ( N log N ) O ( N ) Figure 8: Runtime for solving the Helmholtz equation with heterogeneous wave-speed using GMRESpreconditioned with the the method of polarized traces. The tolerance was set up to 10 − . Left:runtime for solving the low-frequency problem. Right: runtime for solving the high-frequencyproblem with the adaptive basis. . Conclusion In this work we present a numerical method, the ray-FEM, for the high frequencyHelmholtz equation in smooth media based on learning problem specific basis func-tions to represent the wave field. The key information, local ray directions, is ex-tracted from a relative low frequency wave field that has probed the whole domain.These local ray directions are then incorporated into the basis to improve both sta-bility and accuracy in the computation for high frequency wave field. Moreover,both local ray directions and the high frequency wave field can be further improvedthrough more iterations. Numerical tests suggest that our method only requires afixed number of points per wave length without pollution effect as frequency be-comes large. By designing a fast solver for the discretized linear systems an overallcomplexity of order O ( ω d log ω ) is achieved.However, our ray-FEM can not handle singularities of both the amplitude andphase on a mesh. We will develop a hybrid method that combines local asymptoticexpansion near the source and the ray-FEM away from the source in our futurework. 31 cknowledgments Zhao is partially supported by NSF grant (1418422). Qian is partially supportedby NSF grants (1522249 and 1614566).
8. ReferencesReferences [1] I. Babuska, F. Ihlenburg, E. T. Paik, and S. A. Sauter. A generalized finiteelement method for solving the Helmholtz equation in two dimensions withminimal pollution.
Computer Methods in Applied Mechanics and Engineering ,128(3-4):325–359, 1995.[2] I. M. Babuska and S. A. Sauter. Is the pollution effect of the fem avoidablefor the helmholtz equation considering high wave numbers?
SIAM Review ,42(3):451–484, 2000.[3] J.-D. Benamou, F. Collino, and O. Runborg. Numerical microlocal analysis ofharmonic wavefields.
J. Comp. Phys. , 199:714–741, 2004.[4] Jean-David Benamou. An introduction to Eulerian Geometrical Optics (1992-2002).
J. Sci. Comput. , 19(1-3):63–93, 2003.[5] Jean-David Benamou, Francis Collino, and Simon Marmorat. Numerical mi-crolocal analysis revisited.
Research Report, INRIA , 2011.[6] J.-P. B´erenger. A perfectly matched layer for the absorption of electromagneticwaves.
Journal of Computational Physics , 114(2):185–200, 1994.[7] T. Betcke and J. Phillips. Approximation by dominant wave directions in planewave methods. Technical report, 2012.[8] R. Carriere and R. L. Moses. High resolution radar target modeling using amodified Prony estimator.
IEEE Transactions on Antennas and Propagation ,40(1):13–18, Jan 1992.[9] D. Gallistl and P. Peterseim. Stable multiscale Petrov-Galerkin finite elementmethod for high frequency acoustic scattering.
ArXiv e-prints , 2015.[10] E. Giladi. Asymptotically derived boundary elements for the Helmholtz equa-tion in high frequencies.
Journal of Computational and Applied Mathematics ,198(1):52–74, 2007. 3211] E. Giladi and Keller. J. B. A hybrid numerical asymptotic method for scatteringproblems.
Journal of Computational Physics , 174(1):226–247, 2001.[12] R. Hiptmair, A. Moiola, and I. Perugia. A survey of Trefftz methods for theHelmholtz equation.
ArXiv e-prints , 2015.[13] Y. Hua and T. K. Sarkar. Matrix pencil method for estimating parametersof exponentially damped/undamped sinusoids in noise.
IEEE Transactions onAcoustics, Speech, and Signal Processing , 38(5):814–824, May 1990.[14] F. Ihlenburg and I. Babuska. Solution of Helmholtz problems by knowledge-based FEM.
Computer Assisted Mechanics and Engineering Sciences ,(4):397–415, 1997.[15] L.-M. Imbert-Gerard and P. Monk. Numerical simulation of wave propagationin inhomogeneous media using generalized plane waves.
ArXiv e-prints , 2015.[16] H. Jeffreys. On certain approximate solutions of linear differential equations ofthe second order.
Proceedings of the London Mathematical Society , s2-23(1):428–436, 1925.[17] M. Kline and I. W. Kay.
Electromagnetic Theory and Geometrical Optics . In-terscience, New York, 1965.[18] S. Luo, J. Qian, and R. Burridge. Fast Huygens sweeping methods for Helmholtzequations in inhomogeneous media in the high frequency regime.
Journal ofComputational Physics , 270(0):378–401, 2014.[19] J. M. Melenk.
On Generalized Finite Element Methods . PhD thesis, Universityof Maryland, 1995.[20] N. C. Nguyen, J. Peraire, F. Reitich, and B. Cockburn. A phase-based hybridiz-able discontinuous Galerkin method for the numerical solution of the Helmholtzequation.
J. Comput. Physics , 290:318–335, 2015.[21] I. Perugia, P. Pietra, and A. Russo. A plane wave virtual element method forthe Helmholtz problem. 2015.[22] V. F. Pisarenko. The retrieval of harmonics from a covariance function.
Geo-physical Journal International , 33(3):347–366, 1973.3323] L. Rayleigh. On the Propagation of Waves through a Stratified Medium, withSpecial Reference to the Question of Reflection.
Proceedings of the Royal Societyof London Series A , 86:207–226, February 1912.[24] R. Schmidt. Multiple emitter location and signal parameter estimation.
IEEETransactions on Antennas and Propagation , 34(3):276–280, Mar 1986.[25] C. E. Shannon. Communication in the presence of noise.
Proceedings of theIEEE , 86(2):447–457, Feb 1998.[26] C. C. Stolk. A dispersion minimizing scheme for the 3-D Helmholtz equationwith applications in multigrid based solvers.
ArXiv e-prints , 2015.[27] T. Strouboulisa, I. Babuska, and R. Hidajata. The generalized finite elementmethod for Helmholtz equation: Theory, computation, and open problems.
Computer Methods in Applied Mechanics and Engineering , 195(37-40):4711–4731, 2006.[28] H. Wu. Pre-asymptotic error analysis of CIP-FEM and FEM for the Helmholtzequation with high wave number. part I: linear version.
IMA Journal of Nu-merical Analysis , 34:1266–1288, 2014.[29] L. Zepeda-N´u˜nez.
Fast and scalable solvers for the Helmholtz equation . PhDthesis, Massachusetts Institute of Technology, Cambridge MA, USA, 2015.[30] L. Zepeda-N´u˜nez and L. Demanet. The method of polarized traces for the 2DHelmholtz equation.
ArXiv e-prints , 2014.
Appendix A. Stability and error analysis for NMLA
In this section we summarize the stability result and error estimate from [5] forcompleteness. For simplicity we use the single wave case, i.e., N = 1. Moreover, weassume the measurement data is a perturbation to the perfect plane wave data of theform U ( θ ) = U plane ( θ ) + δU ( θ ), where U plane denotes a single plane wave data in theform of (25). Let θ ∗ denote the angle for which θ (cid:55)→ B U ( θ ) is maximum. Assumingthat the noise level satisfies || δU || L ∞ < B ∗ | B | , (A.1)34here B ∗ ≤ .
89 is a pure constant independent of ω and B is the complex amplitudeof the plane wave. Then the error in the angle estimation is given by | θ − θ ∗ | ≤ π L α + 1 ∼ O ( 1 α ) , α = kr ∼ ∞ . (A.2)Similar results can be derived for multiple waves N >
1. We remark that B ∗ (cid:39) .
28, which implies that if the relative noise level does not surpass 28% the angle willbe detected within an error of order O ( kr ). In Benamou’s work [5], an analysis ofa point source shows that | θ − θ ∗ | decreases like O ( ω − / ) when the point x is faraway from the source and the radius of the observation circle is chosen like r ∼ ω − / for large ω . We obtain similar accuracy order for general noisy plane waves undersome smoothness conditions, see details in Appendix B. Appendix B. Error analysis of wave-field as a perturbed plane wave data
As introduced in Section 3.1, NMLA is a tool to process a signal that is (approx-imately) a superposition of plane waves with frequency ω and to extract each planwave component by sampling the signal on a circle/sphere with radius r around areference point. As shown in Appendix A, provided that the perturbation of thesignal is relative small compared to the signal, the estimation of the plane wavedirections converges and the error is O ( ωr ). In this application, we use NMLA toprocess wave-field data, which is the numerical solution to the Helmholtz equation, toextract the directions of dominant wave fronts based on the geometric optics ansatz(7) in the high-frequency regime. Hence it is important to study the wave field dataas a perturbation of plane wave data locally and estimate the error in the ray direc-tions obtained from NMLA. In particular, this analysis allows us to find the optimalchoice of the radius of the sampling circle/sphere, in order to achieve the minimalasymptotic error for the ray direction estimation in terms of the frequency ω of theHelmholtz equation which generates the wave-field data. The result is crucial forboth error analysis and implementation of ray-FEM. Since the wave field data in ourapplication is the numerical solution to the Helmholtz equation, its perturbation canbe composed as the sum of three components:1. numerical error in solving the Helmholtz equation and interpolation error inobtaining data on the sampling circle/sphere for NMLA from the numericalsolution on a fixed mesh,2. the asymptotic error in the geometric optics ansatz,3. the local deviation of a smooth curved wave front from a planar wave front.35n a mesh with mesh size h = O ( ω − ), the last components, which we call thephase error, is the dominant factor among the three. We present below an analysisof the phase error, in which, for simplicity, we only consider one wave front.Let consider a single wave front, u ( x ) = A ( x ) e iωφ ( x ) ; following the notation usedthroughout the paper, assume the reference point to be x , and the small samplingcircle around x to be { x | x − x = r (cid:98) s } , ∇ φ ( x ) = η (cid:98) d , where r (cid:28) , | (cid:98) s | = 1, η = 1 /c ( x ), | (cid:98) d | = 1. A ( x ) = A ( x ) + ∇ A ( x ) · ( x − x ) + O (( x − x ) ) = A ( x ) + r ( ∇ A ( x ) · (cid:98) s ) + O ( r ) ,φ ( x ) = φ ( x ) + ∇ φ ( x ) · ( x − x ) + ( x − x ) T ∇ φ ( x ) ( x − x ) + O (( x − x ) )= φ ( x ) + rη (cid:16)(cid:98) d · (cid:98) s (cid:17) + r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) . Denote φ ( x ) = φ ( x ) + ∇ φ ( x ) · ( x − x ), u ( x ) = A ( x ) e iωφ ( x ) , we have δu ( x ) = u ( x ) − u ( x )= A ( x ) e iωφ ( x ) − A ( x ) e iωφ ( x ) = (cid:2) A ( x ) e iωφ ( x ) + r ( ∇ A ( x ) · (cid:98) s ) e iωφ ( x ) + O ( r ) (cid:3) − A ( x ) e iωφ ( x ) = A ( x ) e iωφ ( x ) (cid:16) e iω [ r ( (cid:98) s T ∇ φ ( x ) (cid:98) s ) + O ( r )] − (cid:17) + r ( ∇ A ( x ) · (cid:98) s ) e iωφ ( x ) + O ( r ) , ∂∂r ( δu ( x )) = ∂∂r (cid:0) A ( x ) e iωφ ( x ) − A ( x ) e iωφ ( x ) (cid:1) = ( ∇ A ( x ) · (cid:98) s + O ( r )) e iωφ ( x ) + A ( x ) e iωφ ( x ) iω (cid:104) η ( (cid:98) d · (cid:98) s ) + r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) (cid:105) − A ( x ) e iωφ ( x ) iωη ( (cid:98) d · (cid:98) s )= ( ∇ A ( x ) · (cid:98) s + O ( r )) e iωφ ( x ) + A ( x ) e iωφ ( x ) iω (cid:104) r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) (cid:105) + (cid:0) A ( x ) e iωφ ( x ) − A ( x ) e iωφ ( x ) (cid:1) iωη ( (cid:98) d · (cid:98) s )= ( ∇ A ( x ) · (cid:98) s + O ( r )) e iωφ ( x ) + A ( x ) e iωφ ( x ) iω (cid:104) r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) (cid:105) + iωη ( (cid:98) d · (cid:98) s ) δu ( x ) . δU ( x ) = (cid:16) iωη ∂∂r + 1 (cid:17) δu ( x )= iωη ( ∇ A ( x ) · (cid:98) s + O ( r )) e iωφ ( x ) + η A ( x ) e iωφ ( x ) (cid:104) r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) (cid:105) +( (cid:98) d · s ) δu ( x ) + δu ( x )= iωη ( ∇ A ( x ) · (cid:98) s + O ( r )) e iωφ ( x ) + η A ( x ) e iωφ ( x ) (cid:104) r (cid:16)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:17) + O ( r ) (cid:105) +( (cid:98) d · (cid:98) s + 1) (cid:110) A ( x ) e iωφ ( x ) (cid:16) e iω [ r ( (cid:98) s T ∇ φ ( x ) (cid:98) s ) + O ( r )] − (cid:17) + r ( ∇ A ( x ) · (cid:98) s ) e iωφ ( x ) + O ( r ) (cid:9) . Hence | δU ( x ) | = (cid:12)(cid:12)(cid:12)(cid:16) iωη ∂∂r + 1 (cid:17) δu ( x ) (cid:12)(cid:12)(cid:12) ≤ |∇ A ( x ) | + O ( r ) ωη + | A ( x ) | η (cid:16) r (cid:12)(cid:12)(cid:12)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:12)(cid:12)(cid:12) + O ( r ) (cid:17) +2 | A ( x ) | ω (cid:16) r (cid:12)(cid:12)(cid:12)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:12)(cid:12)(cid:12) + O ( r ) (cid:17) + 2 r |∇ A ( x ) | + O ( r )= (cid:16) ωη + 2 r (cid:17) |∇ A ( x ) | + (cid:16) | A ( x ) | rη + | A ( x ) | ωr (cid:17) (cid:12)(cid:12)(cid:12)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:12)(cid:12)(cid:12) + | A ( x ) | η O ( r ) + 2 ω | A ( x ) | O ( r ) + O ( r ) . (B.1)As shown in Appendix A, on one hand δU has to be small compared to U . Onthe other hand, the error in direction estimate from NMLA is O ( wr ). Assuming thesmoothness of A ( x ) and φ ( x ), i.e., boundedness of ∇ A ( x ), A ( x ) and ∇ φ ( x ), theleading term in δU is ωr | A ( x ) | (cid:12)(cid:12)(cid:12)(cid:98) s T ∇ φ ( x ) (cid:98) s (cid:12)(cid:12)(cid:12) as ω → ∞ , where (cid:98) s T ∇ φ ( x ) (cid:98) s is thecurvature of the wave front. Hence the radius of the sampling circle can at most bechosen r ∼ O ( √ ω ) as ω → ∞ . Let r = C (cid:15) √ ω , |∇ A ( x ) | ≤ C , | A ( x ) | ≤ C , (cid:12)(cid:12) s T ∇ φ ( x ) s (cid:12)(cid:12) ≤ C , (B.2)Then | δU ( x ) | ≤ C (cid:15) C | A ( x ) | + O (cid:16) √ ω (cid:17) (B.3)37hoose C (cid:15) small enough such that , say 2 C (cid:15) C ≤ , then the perturbation δU ( x )satisfies the condition A.1 for ω large enough, which implies the error in ray directionestimate by NMLA is O ( ω − ). Remark 3.
The above analysis also shows that NMLA can not be used to estimateray directions within a few wavelengths away from the point source since the curvatureof the wave front is of order O ( w ) ..