A dispersion minimizing scheme for the 3-D Helmholtz equation based on ray theory
We develop a new dispersion minimizing compact finite difference scheme for the Helmholtz equation in 2 and 3 dimensions. The scheme is based on a newly developed ray theory for difference equations. A discrete Helmholtz operator and a discrete operator to be applied to the source and the wavefields are constructed. Their coefficients are piecewise polynomial functions of , chosen such that phase and amplitude errors are minimal. The phase errors of the scheme are very small, approximately as small as those of the 2-D quasi-stabilized FEM method and substantially smaller than those of alternatives in 3-D, assuming the same number of gridpoints per wavelength is used. In numerical experiments, accurate solutions are obtained in constant and smoothly varying media using meshes with only five to six points per wavelength and wave propagation over hundreds of wavelengths. When used as a coarse level discretization in a multigrid method the scheme can even be used with downto three points per wavelength. Tests on 3-D examples with up to degrees of freedom show that with a recently developed hybrid solver, the use of coarser meshes can lead to corresponding savings in computation time, resulting in good simulation times compared to the literature.
We consider the discretization on regular meshes of the Helmholtz equation
with large and variable . These methods are widely used for simulations on unbounded domains, for example in exploration geophysics, using domain sizes, in three dimensions, of up to hundreds of wavelengths .
A key issue for such discretizations are the dispersion (phase) errors, that are closely related to pollution errors . Typically, the propagating wave solutions to the discrete and continuous equations have slightly different wavelengths. These wavelength errors are also referred to as phase velocity or phase slowness errors, in which case they are differently normalized. They lead to phase errors in the solution that grow with the distance from the source. A second important consideration is solver cost. The discretized Helmholtz operator should of course be cheap to apply and/or invert.
A class of discretizations, that performs relatively well on these criteria, is given by so called compact finite difference methods, that use a square or cubic stencil in two resp. three dimensions. The corresponding discrete Helmholtz operators can be efficiently applied and inverted compared for example to standard finite difference or finite element methods. Many authors have studied such discretizations and obtained formulae for the coefficients as a function of and the grid spacing . We will discuss these schemes more in detail below, and compare their phase slowness errors with those of standard finite differences and Lagrange finite element methods on regular meshes.
To design such methods, several strategies have been followed. One approach is too construct schemes of higher order, for example order four order six, see  and the references in . Another approach is to stay with second order schemes but minimize the dispersion errors, because these are the dominant errors for long distance wave propagation . From the point of view of phase slowness errors, the sixth order schemes of  and  and the quasi-stabilized FEM (QS-FEM) scheme of  are the best, see the results below. The latter has the smallest phase slowness errors by a substantial margin, but is only available in 2-D.
Alternative methods include higher order finite elements. An advantage of these methods is the better theory for the behavior of the errors in the limit that the grid spacing goes to zero, see for example .
In this paper we introduce a new second order dispersion minimizing scheme in 2 and 3 dimensions with phase slowness errors comparable to those of QS-FEM. Accurate amplitudes are obtained as well using new amplitude correction operators. A theoretical justification is given using a newly developed ray theory for Helmholtz-like difference equations. This theory is remarkably similar to the continuous theory, when both are formulated in terms of the symbols associated with the operators. With numerical examples we show the potential for accurate and fast simulation on relatively coarse meshes. In addition we show applications where the method is used as a coarse level discretization in multigrid solvers
We will briefly describe the methodology and the results. It is known that the second order, compact finite difference discretizations of the Helmholtz operator form a 3 or 5 parameter family, in 2 and 3 dimensions respectively, and that by choosing parameters in a certain way, the phase slowness errors can be reduced compared to standard schemes . When coefficients are allowed to depend on in a piecewise constant  or piecewise linear fashion , they can be further reduced. In this paper we let the parameters depend in a fashion on through third order Hermite interpolation and obtain a further reduction of the phase slowness errors.
Dispersion minimizing schemes are typically intended for use on quite coarse meshes, and a theoretical understanding that does not involve the limit is therefore of considerable interest. For this reason we consider ray theory for Helmholtz-like difference equations.
Ray theory for continuous Helmholtz equations is well known . Solutions are sought in the form . If is smooth, satisfies a certain eikonal equation and a certain transport equation, than such solutions approximate the true solutions increasingly well in the limit . Here we develop a similar theory for Helmholtz-like difference equations. We can then choose the discrete scheme such that the phase and amplitude functions associated with the discrete operator approximate match those of the continuous operator well. As can be expected, schemes with small phase slowness errors have accurate phase functions. By introducing amplitude correction operators, accurate amplitude of the ray-theoretic solutions are obtained.
We are interested in two ways of applying the discretized Helmholtz operators. The first is simply as a discretization of (Equation 1), where the criterion is that the discrete solutions should approximate the true solutions well. Here we are particularly interested in the use of coarse meshes, say downto five or six points per wavelength, which are for example applied in exploration geophysics . The second application is internally in multigrid based solvers. In a multigrid method, the original mesh is coarsened by a factor two one or more times. On each of the new meshes a discretization of the operator is required. In this application the main criterion for a good discretization is that the multigrid method converges rapidly. The results concerning the application in multigrid methods are also of interest for recently developed two-grid or multigrid methods with inexact coarse level inverses , which are currently some of the fastest solvers in the literature. (The method of  will actually be tested here.) Below we will write sometimes the fine level mesh for the original, uncoarsened mesh.
The small phase slowness errors for IOFD suggest that accurate solutions are possible even when quite coarse meshes are used, say downto five or six points per wavelength. We will show that this is indeed the case using numerical examples with constant, and smoothly varying velocity models (recall that with the medium velocity).
We then consider the application of the IOFD discretization as coarse level discretization in multigrid based solvers. We will show that in this case IOFD can be used with very coarse meshes with downto three points per wavelength. With such meshes, solutions are generally not accurate enough for direct use, but the approximate solutions can still be used fruitfully in a multigrid method, where they are refined and iteratively improved. This is established using a set of two-dimensional examples, in which a two-grid method with IOFD at both levels converges rapidly (see also the results discussed in the next paragraph). As explained in , for the good convergence it is necessary to have very small phase slowness errors at these very coarse meshes. The IOFD method (in two and three dimensions) and the QS-FEM method (in two dimensions) are the only discretizations that have this property to our knowledge, and appear to be uniquely suitable for this application.
In 3-D, the fact that a coarser mesh is used does not necessarily imply lower simulation cost. That depends also on the behavior of the solver. To investigate this aspect we present tests with a recently developed solver described in . The solver uses a two-grid method with an inexact coarse level inverse, given by a double sweep domain decomposition preconditioner. As described in the previous paragraph, IOFD will also be used as coarse level discretization. Using the SEG-EAGE Salt Model with up to degrees of freedom as example, we find that for downto six points per wavelength the cost per degree of freedom changes little when the frequency is increased. Computation time compare favorably to some of the results in the literature.
The outline of this work is as follows. In section ? the theory for finite difference discretizations of the Helmholtz equation with constant is developed. The symbols and phase slownesses are defined and the discrete Green’s function is studied. In section ? we consider the case of variable and describe ray theory for discrete Helmholtz equations. In Section 4 we compute the phase slowness errors of various existing schemes, as a reference for the new method. In Section 5 we introduce our new interpolated optimized finite difference method. Section 6 contains some numerical simulations illustrating the accuracy of the solutions when using the IOFD discretization. Section 7 discusses the use of IOFD in multigrid based solvers. Finally, Section 8 contains a brief discussion of some further aspects.
2Theory of discrete Helmholtz equations with constant
In this section we study finite difference discretizations of Helmholtz equation
in case is constant. We will assume the grid is given by . In this and the next section it is convenient to write for multi-indices associated with grid points, such that with is associated the grid point . A difference operator will be viewed as an operator on functions of . In this and the next section the dimension can be any positive integer.
For constant , a finite difference discretization of the Helmholtz operator is a translation invariant difference operator with coefficients depending on the grid spacing and on . By dimensional analysis we may assume that the matrix elements of such a difference operator are defined in terms of a finite set of functions by
where is only nonzero for in some finite set .
We will first consider the action of such an operator in the Fourier domain and define the associated symbol and phase slownesses. We next define a “dimensionally reduced” symbol. Then we consider the discrete Green’s function, i.e. solutions to the equation
where the function on is defined by . We obtain the general solution of this equation in the Fourier domain and the asymptotics in the spatial domain, and determine the same information for the unique outgoing solutions.
In the last part of this section we consider a modification of (Equation 3) where first a function is determined that satisfies
and then is set equal to
In this case we assume and are difference operators of order zero. Based on translation invariance and dimensional reduction, we assume that their matrix elements and are given by
where the , are smooth functions that are only nonzero for in finite sets , . A solution to such a system will be called a modified Green’s function. Our discretization of the Helmholtz equation will be a system of the form (Equation 4) and (Equation 5), where is replaced by the right hand side .
2.1Symbol and phase slownesses
To define the symbol, we first define the forward and inverse Fourier transforms of a function , . They are given by
where the domain of is . For constant the finite difference operator acts like a multiplication in the Fourier domain
where the function , called the symbol, is given by
This is similar to the continuous case, where the Helmholtz operator acts by multiplication with in the Fourier domain.
The Helmholtz equation has propagating plane wave solutions. These are functions that satisfy the homogeneous Helmholtz equation
They are exactly the plane waves for which is in the zeroset of the symbol of (This is of course the set of vectors of length for as defined, but the concept applies more generally.) If is a translation invariant discretization of on we can similarly look for vectors such that
These are the vectors in the zero set of .
If is a discretization of the Helmholtz operator then typically the set is close to, but not identical to . In other words, there are small differences in the wave vectors of the propagating waves, . If and can be parameterized by angle, i.e.
and similar for (for this is of course the case), we define the relative wave number error as a function of by
Closely related quantities are the phase slowness and phase velocity errors. If , , there are associated phase slowness and phase velocity vectors given by
The actual phase error between a numerical and an exact solution is given by (see also subsection ?)
where is the distance between source and observation point, is the wavelength and is the phase slowness error associated with the particular angle of propagation. Because it is proportional to the phase error easily may become dominant if is not careful in the choice of discretization in the high-frequency regime.
2.2Dimensional reduction and Helmholtz-like symbols
In case of coefficients, the symbol for arbitrary can be expressed in terms of that for
By dimensional reduction the symbol , can be written as
and similar for .
To obtain the results below, we assume that the symbol is Helmholtz-like as defined in the following
It follows that is Helmholtz-like if is Helmholtz like.
2.3The discrete Green’s function
A Green’s function for the discrete equation will be defined as a solution of (Equation 3). If is constant the equivalent equation for the Fourier transform reads
We will first describe the general solution to this equation in the Fourier domain. Then we will consider the asymptotics in the spatial domain. Using the results obtained, we can then derive a unique outgoing Green’s function in the Fourier domain, and state its asymptotics.
Due to the zeros of , problem (Equation 12) has non-unique, distributional solutions. To explain their nature, we recall the closely related one dimensional problem to determine all such that
where is a free constant. Here the distribution is defined by
In other words is in the kernel of , while is a particular solution to (Equation 13).
In case of (Equation 12) we similarly have a nonzero kernel of , with elements , where denotes the singular function of , which is the distribution given by
and is any distribution on . For functions on with zero set such that for all , the principal value can be defined as follows. Let , then
The freedom in the choice of is related to the fact that in the spatial domain one can add any linear combination of plane waves with and still have a solution.
Let be the inverse Fourier transform of a solution for some smooth
We will study the asymptotic behavior of this integral for large using the method of stationary phase. For , we define a certain curvature-like quantity as follows. After rotating the coordinates, we may assume that is parallel to the -th coordinate axis and that is locally a graph
By the assumptions has a nondegenerate local maximum at . Let denote the eigenvalues of the second derivative matrix . We define
For this is the Gaussian curvature of the surface. In the following proposition denotes the inverse function of the map defined in ( ?). The result and its proof have some similarity with results of Lighthill .
We start with the first integral in (Equation 15). For the domain is really a torus and the integrand is as a function on the torus. It is convenient to replace the integral on the torus by an integral over a bounded subset of . Let be a smooth, positive function supported in , that is one on and satisfies for , and let
Then we can write
for , where is the periodic extension of , and this formula may also be considered for . We assume is sufficiently small such that is supported in .
We will write , and consider the limit . We assume coordinates are rotated such that , using the same notation for the new coordinates as used so far for the old coordinates.
The integral on the right hand side of (Equation 18) will be written as a sum of integrals over subsets using a partition of unity. For some smooth cutoff function , denote
We may assume there are four different types of
is one on a neighborhood of
is one on a neighborhood of
on we can write as a graph
We consider these four cases in the limit using the method of stationary phase . In case (iv) the integral for any by the lemma of non-stationary phase and we don’t need to consider this case further. In case (iii) we can write
for some smooth function and perform the integral over . This yields a smooth function of . By the lemma of non-stationary phase it follows that again .
In case (i) we can write locally as a graph . For brevity denote We observe that we can write
where are smooth, compactly supported functions and around and . Then for the first term the lemma of non-stationary phase can be invoked. Hence this term is for any . For the third term, the same result can be obtained using integration by parts. For the second part we recall the standard Fourier transform , it follows that
As a consequence, we obtain
any . We can now apply the stationary phase lemma. The function has its maximum at and can be expanded as, possibly after a further rotation of coordinates
see the discussion preceding (Equation 17). This yields
The contribution in case (ii) can be computed similarly, resulting in
For the surface integral
we again assume and consider the limit . A partition of unity is applied and, by the method of stationary phase, the only contributions that are not for any come from neighborhoods of . To determine the contribution from a neighborhood of , we assume that is parallel to the -th coordinate axis so that is locally given by a graph , . The method of stationary phase can be applied directly. The only contributions come from neighborhoods of , and can be computed similarly as for the integral (Equation 19).
The contribution , and the two contributions from the integral (Equation 20) together give the result.
We state this as a theorem and include the asymptotic expression for the solution in the result.
The above analysis can be repeated for continuous, Helmholtz like operators with the same result (see also ). For the usual Helmholtz operator in dimensions we have that given by , , and
so that the highest order asymptotic expansion actually equals the well known outgoing Green’s function.
2.4The modified discrete Green’s function
It follows from theorem ? that if and have the same zero sets, i.e. identical phase slownesses, then the solutions to and have asymptotically the same phase. The amplitudes however will differ by a factor evaluated at the zero set. In this subsection we consider therefore the solutions to the equations (Equation 4) and (Equation 5), which, as we will see, obtain different amplitudes.
The Fourier transformed solution to (Equation 4) and (Equation 5) is given by the product of the solution given in proposition ? and a factor . Using this, we can formulate a result similar to Theorem ?. In this case the adjective outgoing refers to the solution of (Equation 4). The result can be proven by similar arguments as used to prove proposition ? and theorem ?.
Summarizing our findings so far, the discrete solutions to (Equation 4) and (Equation 5) are asymptotically equal to the solutions of the continuous Helmholtz equation if the following two conditions are satisfied
and have the same zero sets
for all such that
3Theory of discrete Helmholtz equations with variable
In this section we define a class of discrete approximations to the Helmholtz operator with variable , together with the associated symbols. This is the topic of subsection ?. We then study ray-theoretic solutions to the equation (Equation 3) and to the set of equations (Equation 4), (Equation 5), where and are now variable coefficient operators.
We assume that , where is smooth and we consider the limit . In the discrete case we assume that . Ray-theoretic solutions are then based on the ansatz
for some smoothly varying and . For the continuous Helmholtz equation, such solutions are well known and are constructed in two steps. First the ansatz (Equation 22) is inserted in the PDE, and an expansion in is performed. Requiring that the highest order terms vanish leads to the eikonal equation for and the transport equation for . Secondly, initial/boundary conditions for these equations are obtained from the asymptotic behavior of the constant coefficient solutions. In this way, the solution modulo an error of lower order in is obtained.
For our class of difference equations we follow the same program. The constant coefficient solutions were already analyzed in subsection ?. In subsection ? we find a nonlinear first order PDE for and a transport equation for . Remarkably, we obtain the same equations in the continuous and discrete case when formulated in terms of the symbols (which are defined for both continuous and discrete problems). See  and  for the continuous case and methods used in that case as well as here.
In the last part of this section we consider the ray-theoretic solutions to (Equation 4) and (Equation 5). The conditions ( ?) and ( ?) from subsection ? for and to obtain accurate solutions, need to be modified and extended to have the same ray-theoretic phase and amplitude in the continuous and discrete case. The operator should be discretized using a symmetric discretization (with , see below) and we should have . This is the topic of subsection ?.
3.1Symbols and operators for variable
In case depends on , finite difference discretizations of the Helmholtz operator may depend in different ways on the function . For example, the coefficients may depend on and its derivatives at , but they may also depend on at different points, for example on and . We will consider a class of difference operators , where the matrix elements depend only on the value of at , where is a fixed constant. In other words we consider operators with matrix elements of the form
Note that the operator is symmetric if and . This will turn out to be an appropriate choice for a discrete Helmholtz operator. We will assume that is defined for all , not only those in the grid. Similar we assume that for we have
and similar for .
For such operators it is not obvious how to define the symbol. To find an appropriate definition, we first consider how to define an operator from a symbol in the continuous case. This is the subject of pseudodifferential operator theory, and can be done with the formula 
A map from a function to an operator such as is called a quantization. For , the previous formula is the standard left-quantization, is the right quantization and is the Weyl quantization. If , then , independently of which of these quantizations is used.
To obtain a symbol associated with the operator defined in (Equation 23) we rewrite the expression for as follows
Using the Fourier domain representation this can be rewritten as
This can be written in similar form as (Equation 24), namely as
With these definitions, the symbol (Equation 26) for variable coefficients may also be expressed entirely in terms of
A symbol is called Helmholtz like if it satisfies the definition for each fixed .
3.2Ray-theoretic equations for amplitude and phase
In this section we consider the high-frequency limit . We assume that , and recall that , where is . The operator and the symbol become -dependent. By we denote the symbol for .
For other values of we find that
We consider the action of on functions of the form
where and are functions. From the symbol and the phase function one can derive naturally a vector field, which we call (cf. )
This vector field is determined by as follows
The proof uses a Taylor expansion of the phase function to second order
a Taylor expansion of the amplitude to first order
and a Taylor expansion of the matrix coefficients to first order
The exponent is then written as a product of three factors
These expansions are inserted in the sum
The factor can be put in front of the expression outside the summation
We next use the expression for the symbol , the following expressions for the derivatives of
and an expression for the derivatives of
The result is similar to the result in Proposition 4.3.2 of . To find and such that
the phase function must satisfy the equation
which is a nonlinear first order equation like an eikonal equation, and the amplitude must satisfy a transport equation
For , this equation conserves .
Ray theoretic solutions to equation (Equation 3) can now be constructed just as in the continuous case. By a rescaling, theorem ? can be used to obtain the asymptotics of a solution for and . The amplitude and phase from formula ( ?) can hence be used as initial/boundary values for the eikonal equation for and the transport equation for , and these and can be determined from these equations, where we note that the eikonal equation may not have globally defined solutions, just as in the continuous case.
We briefly recall the continuous equivalent of Proposition ?. The following is basically a reformulation of proposition 4.3.2 of  and can be proven using the method of stationary phase found in the same text.
In this section we consider ray-theoretic approximations and for the solutions and to (Equation 4) and (Equation 5). Assume we have reference ray-theoretic solutions associated with Helmholtz equation where , i.e. satisfies the eikonal equation and the transport equation with appropriate initial conditions. In the following we let , .
and have the same zero sets
and are identical and satisfies
for all such that ;
and are derived from their respective symbols using quantization .
We argue that in this case, to highest order has the same ray-theoretic approximation as . We omit a formal proof, because the arguments are similar as those used above.
The construction of the phase and amplitude functions and proceeds almost in the same way as for solutions to (Equation 3). Eikonal and transport equations are as follows from Proposition ?. The constant coefficient solutions differ by a factor and have the same phase, resulting in different initial/boundary conditions, such that on a small sphere around 0, where we impose the initial/boundary conditions for the eikonal and transport equations, we have
As a result, everywhere. While we have different transport equations the operators and are scaled versions of each other
It follows from this fact and the transport equation for , that
This and (Equation 29) shows that
The function is given by applying to . The action of on is to highest order equal to a multiplication by , so that
concluding the argument.
4Phase slowness errors for existing discretizations
In this section we will describe three types of discretizations of the Helmholtz equation (Equation 1), namely standard finite differences, compact finite differences and Lagrange finite elements on regular meshes. We then compute phase slowness errors to compare the performance of the different methods in this respect, and to obtain reference values for our new method constructed below.
We modify the notation compared to the previous two section. In this section the degrees of freedom for all three types of methods are denoted by (in three dimensions) and associated with a regular mesh with grid spacing . For finite element methods of order the cells are of size , and cell boundaries are located at . Occasionally we will use or to denote the dimension of space.
4.1Standard finite differences
In a standard finite difference discretization of the operator each of the one-dimensional second derivatives in the Laplacian is approximated by a central difference approximation of the given order. These are given by
where the are as in the following table for