# Preconditioned Krylov Subspace Methods for Sixth Order Compact Approximations of the Helmholtz Equation

Authors: Yury Gryazin
Journal: ISRN Computational Mathematics
 (2014)

---

## Abstract

We consider an efficient iterative approach to the solution of the discrete Helmholtz equation with Dirichlet, Neumann, and Sommerfeld-like boundary conditions based on a compact sixth order approximation scheme and lower order preconditioned Krylov subspace methodology. The resulting systems of finite-difference equations are solved by different preconditioned Krylov subspace-based methods. In the analysis of the lower order preconditioning developed here, we introduce the term “kth order preconditioned matrix” in addition to the commonly used “an optimal preconditioner.” The necessity of the new criterion is justified by the fact that the condition number of the preconditioned matrix in some of our test problems improves with the decrease of the grid step size. In a simple 1D case, we are able to prove this analytically. This new parameter could serve as a guide in the construction of new preconditioners. The lower order direct preconditioner used in our algorithms is based on a combination of the separation of variables technique and fast Fourier transform (FFT) type methods. The resulting numerical methods allow efficient implementation on parallel computers. Numerical results confirm the high efficiency of the proposed iterative approach.


---

## Body

## 1. Introduction

In recent years, the problem of increasing the resolution of existing numerical solvers has become an urgent task in many areas of science and engineering. Most of the existing efficient solvers for structured matrices were developed for lower-order approximations of partial differential equations. The need for improved accuracy of underlying algorithms leads to modified discretized systems and as a result to the modification of the numerical solvers (see, e.g., [1]).The use of a lower-order preconditioner for efficient implementation of high-order finite-difference and finite-element schemes has been under consideration for a long time (see, e.g., [2, 3]). In this paper, a compact sixth order approximation finite-difference scheme is developed, and a lower-order approximation direct solver as a preconditioner for an efficient implementation of this compact scheme for the Helmholtz equation in the Krylov subspace method framework is considered. This approach allows us to utilize the existing lower-order approximation solvers which significantly simplifies the implementation process of the higher resolution numerical methods.The model problem considered in the paper is the numerical solution of the Helmholtz equation(1)∇2u+k2u=f,inΩ,
with the Dirichlet, Neumann, and/or Sommerfeld-like (radiation) boundary conditions
(2)u=0,on⁡∂Ω1,un=0,on⁡∂Ω2,un-iku=0,on⁡∂Ω3,
where Ω={0≤x,y,z≤a}, k is a complex valued constant, and ∂Ω1, ∂Ω2, and ∂Ω3 are different boundary sides of the rectangular computational domain Ω.It is known that for a given error level in the numerical approximation to the solution of the Helmholtz equation, the quantity(Re(k))p+1hp needs to be constant, where p is the order of finite-difference scheme and h is the grid size. This phenomenon is known as “pollution” [4, 5]. One way of reducing the pollution error is to increase the order of accuracy of the scheme. In this paper a sixth order compact finite-difference scheme is considered to address this problem. In the cases of the 3D Helmholtz and Dirichlet boundary conditions, this scheme was proposed by Sutmann in [6]. The 2D version of this scheme with Dirichlet and Neumann boundary conditions was developed in [7, 8]. In this paper, the method is extended to the explicit compact sixth order approximation of Neumann and Sommerfeld-like boundary conditions in the 3D case. The extension of the known approach for the Neumann boundary conditions is straightforward. But in the case of the Sommerfeld-like boundary conditions the method is nontrivial and requires the introduction of a new auxiliary function. In the case of the variable coefficient Helmholtz equation and Sommerfeld-type boundary conditions, the fourth order compact approximation scheme was considered in [9].The resulting discretization leads to a system of linear equations with block 27-diagonal structure. In general, the matrix of this system is neither positive definite nor Hermitian. Hence, most iterative methods either fail to converge or converge too slowly, which is impractical. For the solution of this system we propose using a combination of Krylov subspace-based methods and the FFT preconditioner. Concerning other approaches for the solution of the problem (1)-(2), we refer to Bayliss et al. [10] for a preconditioned conjugate-gradient algorithm, Kim [11] for a domain decomposition method, and Douglas et al. [12] for an ADI algorithm. A very efficient multigrid method based on a shifted-Laplace preconditioner for the Helmholtz equation is presented in [13]. The analysis of the multigrid algorithms in the case of nonsymmetric and indefinite problems can be found in [14].On the other hand, the solution of this problem by a direct method based on Gaussian elimination requires a prohibitive amount of additional storage and computer time and thus has limited use. The most promising results in the solution of a similar problem have been obtained by preconditioned Krylov subspace methods in [15]. In this paper we generalize some approaches developed for the second-order central difference discretization of the Helmholtz equation by Elman and O’Leary [16, 17], the author, and others [18] to the case of compact sixth order approximation scheme.The key to the fast convergence of the suggested iterative method is the choice of the preconditioning matrix. In our algorithm we use a preconditioner based on the second order central difference approximation of the Helmholtz equation and corresponding boundary conditions. The inversion of the preconditioning matrix at each step of the Krylov subspace method is done by a direct solver based on the FFT technique (see, e.g., [19]) which requires O(N3log⁡N) operations, where N is the number of grid points in each direction.Numerical experiments with test problems demonstrate the high resolution of the sixth order compact finite-difference scheme, as well as the computational efficiency of our preconditioned Krylov subspace-type numerical method. In most situations, the GMRES method demonstrates the best convergence properties. However, in the case of sufficiently fine grids, we propose using a simple stationary two-level method (see, e.g., [20]). This method was naturally constructed in the analysis of the convergence of the GMRES algorithm. So, the choice of the parameters in this method is not based on the minimization of the spectral radius of the iteration matrix on each step but rather on the construction of a particular linear combination from a Krylov subspace. This is why, for the purpose of this paper, we took the liberty of calling this approach the simplified Krylov subspace (SKS) method. On sufficiently fine grids this method requires much less processor time than the GMRES method, though the number of iterations until convergence is still larger than in the case of the GMRES method as should be expected. Also, since the implementation of SKS algorithm does not require the calculation of scalar products, it has greater potential than the GMRES method for implementation on parallel computers. A distinguishing feature of these approaches is that the number of iterations required for the convergence of Krylov subspace iterations decreases as the size of the discretized system increases. To explain this fact we introduce “the order of the preconditioned matrix” as a parameter to quantify the rate of convergence of the condition number of the preconditioned matrix to 1 on a sequence of grids. In some simple situation, we were able to find this parameter analytically. We believe that this parameter is more informative than the commonly used ‘‘an optimal preconditioner” (see, e.g., [21, page 196]) in the case of lower-order preconditioners and may be used as guide in the further development of preconditioners of similar type. But we must notice that even in this paper this parameter has limited application.This conclusion is based on the theoretical analysis of some simple situations and confirmed in our numerical experiments. Some preliminary numerical results on the application of this approach were presented at the 10th International Conference on Mathematical and Numerical Aspects of Waves, Vancouver, Canada, 24–29 July, 2011 [22].The rest of the paper is organized as follows. In Section2, the main idea of the proposed sixth order approximation compact method is presented in the case of the 1D problem. To analyze the convergence of the algorithms developed here, variations of known convergence estimates for the Krylov subspace methods are considered. In this section, simplified approaches based on the Krylov subspace iterations are also presented. Section 3 focuses on the development of the sixth order compact approximation scheme in the case of Neumann and Sommerfeld-like boundary conditions. In Section 4, the effectiveness of the proposed algorithms is demonstrated on a series of test problems.


## 2. A One-Dimensional Model Problem

LetA and Ap be matrices derived from six and second order approximations to (1) and the Dirichlet boundary condition (2) using a mesh xi=ih,i=0,1,…,N+1, and h=1/(N+1). The standard notation for the first and second order central differences at ith grid point is given by
(3)δxui=ui+1-ui-12h,δx2ui=ui+1-2ui+ui-1h2,
where ui=u(xi). The finite difference operators δy, δz, δy2, and δz2 used in the following sections are defined similarly. The sixth order approximation to the second derivative at the ith grid point can be written as
(4)ui′′=δx2ui-h212ui(4)-h4360ui(6)+O(h6).
As usual, in the case of compact schemes, we need to use the original equation to find appropriate relations to eliminate the fourth and sixth derivatives in (4). By using the second order central difference of the fourth derivative of u at x=xi(5)δx2ui(4)=ui(6)+O(h6)
and the expression for the fourth derivative of u(x) from (1)
(6)ui(4)=-k2ui′′+fi′′,
the relation (4) can be expressed in the form
(7)(1-k2h212(1+h230δx2))ui′′=δx2ui-h212(1+h230δx2)fi′′+O(h6),
where fi=f(xi). After using (7), the discretized system corresponding to the compact sixth order approximation scheme for (1)-(2) can be presented as
(8)d1Ui-1+d2Ui+d3Ui+1=Fi,i=1,…,N,U0=0,UN+1=0,Fi=h2(1-7k2h290)fi-k2h4360(fi-1+fi+1)+7h490fi′′+h4360(fi-1′′+fi+1′′),d1=d3=1-k4h4360,d2=-2+k2h2-7k4h490,
where Ui is the sixth order discrete approximation to u(xi). This system can also be rewritten in the form AU=F, where
(9)A=d1Λ+d2I,(10)Λ=[01⋯0101⋮⋱⋱⋱⋮1010⋯10].
In a similar way the preconditioning matrix based on the second order central difference approximation can be presented as
(11)Ap=Λ-(2-k2h2)I.
Finally the right preconditioned system can be written as
(12)AAp-1Y=F,ApU=Y.The main goal of this paper is to demonstrate the computational efficiency of a preconditioning technique based on the use of a lower-order approximation discrete system in the numerical implementation of higher order compact finite-difference schemes. To consider this construction in a more general form, we introduce the definition of akth order preconditioned system.Definition 1.
Given two nonsingularN×N matrices A and Ap, a system in the form of (12) is said to be a kth order preconditioned system if the N×N matrix AAp-1 is diagonalizable and can be expressed as V-1(I+hkD)V, where h<1; D=diag⁡(d11,…,dnn) with max⁡i|dii|<M, where M does not depend on h and V is an N×N matrix of eigenvectors of AAp-1.Using this definition we present the convergence analysis of the GMRES method applied to the system (12) in the following theorem.Theorem 2.
Let a system in the form of (12) be a kth order preconditioned system. Then the nth iteration U(n) of the GMRES method applied to this system satisfies the convergence estimate
(13)∥rn∥2≤κ2(V)(Mhk)n∥r0∥2,
where κ2(V)=∥V-1∥2∥V∥2 and rn=F-AU(n).Proof.
Since the matrixAAp-1 is diagonalizable, it is well known (see, e.g., [23, 24]) that the residual of nth iteration of GMRES algorithm satisfies
(14)∥rn∥2≤κ2(V)min⁡p∈Pn,p(0)=1max⁡i=1,…,N|p(1+hkdii)|·∥r0∥2.
Consider the polynomial of nth degree p^n(x)=1-xyn-1(x), where yn-1(x)=∑k=0n-1‍αk+1nxj and the coefficients αkn satisfy p^n(1+x)=-αnnxn,n≥1. It is easy to see that the coefficients αkn,k=1,…,n satisfy
(15)∑k=1n‍αkn=1,∑k=l-1n‍(kl-1)αkn=0,l=2,…,n,
and the unique solution of this system is given by αkn=(-1)k-1(nk),k=1,…,n. Then the convergence estimation for the GMRES method becomes
(16)∥rn∥2≤κ2(V)min⁡p∈Pn,p(0)=1max⁡i=1,…,N|p(1+hkdii)|·∥r0∥2≤κ2(V)max⁡i=1,…,N|p^n(1+hkdii)|·∥r0∥2≤κ2(V)(Mhk)n∥r0∥2.The proof is somehow trivial and could be significantly simplified but we present it in such a form to use later in the construction of the simplified iterative method based on this proof.In the same way we can derive another useful estimate expressed in the following corollary.Corollary 3.
Let the matrixAAp-1 in the system (12) be expressible as V-1(I+D)V, where D=diag⁡(d11,…,dnn), M=max⁡i|dii|<1, and V is the N×N matrix of eigenvectors of AAp-1. Then the nth iteration Un of the preconditioned GMRES method applied to this system satisfies the convergence estimate
(17)∥rn∥2≤κ2(V)Mn∥r0∥2.We can derive another well-known estimate based on a standard application of Chebyshev polynomials (see, e.g., [24]). We consider only the case of real valued matrices. The following theorem provides the details of the proposed technique.Theorem 4.
Let the matrixAAp-1 in the system (12) be diagonalizable and expressible in the form V-1(I+D)V, where D=diag⁡(d11,…,dnn), -1<m^=min⁡idii≤M^=max⁡idii, and V is the N×N matrix of eigenvectors of AAp-1. Then the nth iteration Un of the preconditioned GMRES method applied to this system satisfies the convergence estimate
(18)∥rn∥2≤2κ2(V)(M^-m^4(1+m^))n∥r0∥2.Proof.
This estimation immediately follows from the general complex valued matrix result (see, e.g., [24, page 206]).Next we will apply the developed convergence estimations to the analysis of the solution of the system (12).Theorem 5.
Givenk, let N×N matrices A and Ap be defined by (9) and (11) with h<2π/10k and min⁡j=1,…,N|(4sin2(jπh/2)/h2)-k2|≥k2δ0>0, where δ0 is independent of h. Then AAp-1 is a second order preconditioned system and one has the convergence estimate
(19)∥rn∥2≤(k2h212δ0)n∥r0∥2.Proof.
The matricesA and Ap have the same set {Vj,j=1,…,N} of eigenvectors as the matrix Λ from (10). The eigenvectors are given by (see, e.g., [17])
(20)Vjl=2hsin(jlπh),j,l=1,…,N.
The matrixV is unitary, so ∥V-1∥2=∥V∥2=1. The eigenvalues of the matrices Ap and A are given by
(21)λjp=-4sin2(jπh2)+h2k2,forj=1,…,N,λj=-4sin2(jπh2)+h2k2-h4k4(14+cos⁡⁡(jπh))180,λj=-4sin2(jπh2)+h2k20000000-forj=1,…,N,
respectively. We can write the eigenvalues of the preconditioned matrix AAp-1 in the form
(22)λjλjp=1+h4k4(14+cos⁡(jπh))/1804sin2(jπh/2)-h2k2,forj=1,…,N.
The second term in the right-hand side can be estimated using m≤|(λj/λjp)-1|≤M, where m=13h4k4/720 and M=h2k4/(12δ0). So AAp-1 is a second order preconditioned matrix and using (13) we obtain the estimate stated in the theorem.The first condition of this theoremhk<2π/10 is a common natural restriction that comes to play even when one wants to visualize a numerical solution for qualitative analysis. It is just not an accurate representation of the solution if there are fewer than 5 points per half wavelength. Moreover, to avoid the so-called “pollution” phenomenon Bayliss et al. [5] introduced the restriction that kp+1hp should be constant; that is, to maintain a fixed accuracy of a scheme, the number of grid points must grow as k1/p, where p is the accuracy order of the scheme. So, if one satisfies the restriction that avoids “pollution” phenomenon, the first condition of the theorem is satisfied almost automatically and it is not much of a restriction at all.The second condition of the theorem is a requirement that avoids the discrete spectrum of the operator. Since the resolvent set of a bounded linear operator is open, we can always do this by introducing a small change tok2 if necessary.Next, we present how we can chooseδ0 in some important cases.In the simplest case ofk2<9, we can take δ0=1/27.Consider how to get an estimate forδ0 in the case of a given k>π and a sequence of grids with the grid sizes h1>h2>h3⋯. The first restriction of the theorem gives us the condition kh<2π/10. We can see that if the j in the argument of min⁡j|(4sin2(jπh/2)/h2)-k2| was to take on values of all real numbers from 1 to N, then the minimum of the expression under consideration would be zero.Now let us denoteα(j)=jπh/2 and consider for which α the expression is zero if kh=2π/10. We consider the worst-case scenario to justify this approach for all kh<2π/10. It is trivial to find that α0=sin-1(kh/2)=sin-1(π/10)≈.32. So, it is clear that the minimum of the function occurs either at j0=⌊2α0/πh⌋ or at j1=j0+1 since sin(α(j)) is increasing on the interval 1<j<N. There exists a neighborhood of α0 which includes both α(j0) and α(j1) such that α<.64 which allows us to use the Taylor series to represent sin2(α) in this neighborhood. We use only the first two terms in the series; that is, sin2(α)=α2-(1/3)cos⁡(2ξ)α4, 0<ξ<0.64. Now we can substitute this expression
(23)min⁡j=1,…,N|4sin2(jπh/2)h2-k2|=min⁡j=j0,j1|j2π2-112cos⁡(2ξ)(jπ)4h2-k2|≥k2min⁡[1-j02π2k2,j12π2k2-1-112(j1π)4(h1k)2]≥k2δ0.
So, we just need to avoid the values for k and h1 for which one of the terms in the brackets is nonpositive. In some situations we would need to choose a sufficiently small grid size step h1 for the coarsest grid in a sequence. For example, in the cases k=20 and h1=1/32, min⁡(1-62π2/400,72π2/400-1-(1/12)((7π)4/(400·322)))>min⁡(.11,.1575)=δ0. This δ0 is independent of h and can be used in our convergence analysis on a sequence of grids.It follows from (19) that the number of iterations for the Krylov subspace algorithm decreases as the grid size h decreases. This result proves that the algorithm developed here yields a very effective iterative technique for the implementation of a higher order approximation scheme. Later we will consider the extension of this method to 3D problems.The proof of Theorem2 suggests a simplified version of the Krylov subspace method for the solution of (12). Indeed, using the expression for the exact solution of the system (15), we have αnn=-αn+1n+1, for n=1,… and the following recurrence expression for yn(24)yn(z)=1-zyn-1(z)+yn-1(z),wheren≥1.
This allows us to derive the SKS algorithm described in the table labelled Algorithm 1.Algorithm 1:Simplified Krylov subspace (SKS) algorithm.
(1) LetY(0)=Ap-1U(0)=0 be the initial approximation, let the tolerance be tol,
and let the maximum number of iterations beM
(2)r0=F, Y=r0 and δ0=∥r0∥2
(3)j=1
(4)while  j<M and ε≥tol  do
(5)U(j)=Ap-1Y
(6)w=AU(j)
(7)r=r0-w
(8)δ=∥r∥2
(9)ε=δδ0
(10)Y=Y+r
(11)j=j+1
(12)end while
(13)U(j) is the iterative solution of AU=F.This is a simple stationary two-level method (see, e.g., [20]). However, the choice of the parameters in this algorithm is not based on the minimization of the spectral radius of the iteration matrix on each step but rather on the construction of a particular linear combination of vectors from a Krylov subspace. This is why, for the purpose of this paper, we call this approach the simplified Krylov subspace (SKS) method. The proposed method does not require an orthogonalization procedure; that is, it requires no inner products, which is attractive for an implementation in a parallel computing environment. In addition, this method does not use estimates of maximum and minimum eigenvalues of the preconditioned matrix which is the drawback of the Chebyshev acceleration algorithm (see, e.g., [24]). The Chebyshev acceleration method is naturally constructed in the proof of Theorem 4 and we will compare the numerical effectiveness of the methods in Section 4. We should notice that in most numerical experiments, the GMRES method exhibits the best convergence properties, but in some situations the SKS algorithm proves to be more efficient. Some such examples are considered in Section 4.The SKS algorithm also gives a good criterion for evaluating the quality of a preconditioner in the form of the order of preconditioned matrix (12). In this paper we will calculate this as follows. Consider two grids with the grid sizes h and γh, where γ>1. Let the l2-norm of the residuals of the SKS method on the first two iterations be ∥r1h∥2 and ∥r2h∥2 on the first grid and ∥r1γh∥2 and ∥r2γh∥2 on the second grid. Let εh=∥r2h∥2/∥r1h∥2 and εγh=∥r2γh∥2/∥r1γh∥2. Now we can approximately calculate the order of preconditioned matrix by
(25)ψ=ln⁡εγh-ln⁡εhln⁡γ.
We consider this parameter in the discussion of the results of the numerical experiments.


## 3. Three-Dimensional Problems

### 3.1. A Sixth Order Approximation Compact Scheme

In this section we present a three-dimensional compact sixth order approximation finite-difference scheme that was first introduced in [6], and we develop a sixth order compact explicit approximation of the Neumann and Sommerfeld-like boundary conditions (2). In this discussion we consider a uniform grid Ωh={(xi,yj,zk)∣xi=ih,yj=jh,z=kh;i,j,k=0,…,N-1,h=a/(N-1)}.First, we consider a finite-difference approximation of (1) in the form
(26)δx2ui,j,k+δy2ui,j,k+δz2ui,j,k+k2ui,j,k+Ti,j,k=fi,j,k,
where ui,j,k=u(xi,yj,zk) and fi,j,k=f(xi,yj,zk) and
(27)Ti,j,k=-[h212(∂4u∂x4+∂4u∂y4+∂4u∂z4)hhh+h4360(∂6u∂x6+∂6u∂y6+∂6u∂z6)]i,j,k+O(h6).Using the appropriate derivatives of (1) we can write the sixth order compact approximation of the Helmholtz equation in the form
(28)(δx2+δy2+δz2)Ui,j,k+h26(1+k2h230)×(δx2δy2+δx2δz2+δy2δz2)Ui,j,k+h430δx2δy2δz2Ui,j,k+k2(1-k2h212+k4h4360)Ui,j,k=(1-k2h212+k4h4360)fi,j,k+h212(1-k2h230)∇2fi,j,k+h4360∇4fi,j,k+h490(∂4f∂x2∂y2+∂4f∂x2∂z2+∂4f∂y2∂z2)i,j,k.
For the two-dimensional case this scheme was proposed in [7]. In the three-dimensional case with Dirichlet boundary conditions, it was developed in [6]. We will consider the iterative implementation of this scheme based on preconditioned Krylov subspace-type algorithms.


### 3.2. Boundary Conditions

The implementation of the Dirichlet boundary conditions (2) is straightforward but the explicit compact approximation of Neumann and Sommerfeld-like boundary conditions requires careful consideration. The sixth order compact approximation of the Neumann boundary condition in the two-dimensional case was considered in [7]. Here we extend this approach to the three-dimensional problem.Neumann Boundary Conditions. In this subsection the Neumann boundary conditions are considered in the form
(29)∂u∂z|z=0=β(x,y).We restrict our consideration to only one side of the computational domainΩ. By using Taylor series, we can derive the sixth order approximation formula
(30)δzui,j,0=∂u∂z|i,j,0+h26∂3u∂z3|i,j,0+h4120∂5u∂z5|i,j,0+O(h6).
To express the third and fifth derivatives in (30) we can differentiate the original Helmholtz equation (1), assuming sufficient smoothness of the solution and the right-hand side. After substitution of these derivatives into (30), the resulting expression is
(31)δzui,j,0=[∂u∂z]i,j,0+h26×[∂f∂z-∂3u∂x2∂z-∂3u∂y2∂z-k2∂u∂z]i,j,0+h4120[∂3f∂z3-∂5u∂z3∂x2-∂5u∂z3∂y2hhhhhhh-k2(∂f∂z-∂3u∂x2∂z-∂3u∂y2∂z-k2∂u∂z)]i,j,0+O(h6).Next, we approximate the third order mixed derivatives in the second bracket of (31) by a fourth order approximation formula and in the second line using a second order approximation scheme, which yields
(32)δzui,j,0=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0-h26[δzδx2u+δzδy2u-h26(∂5u∂x2∂z3+∂5u∂y2∂z3)hhhhhh-h212(∂5u∂x4∂z+∂5u∂y4∂z)]i,j,0+h4k2120[δzδx2ui,j,0+δzδy2ui,j,0]-h4120[∂5u∂x2∂z3+∂5u∂y2∂z3]i,j,0+O(h6).
Now we can simplify this and write
(33)δzui,j,0+h26(1-h2k220)[δzδx2ui,j,0+δzδy2ui,j,0]=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0+7h4360[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0-7h4180[∂5u∂x2∂y2∂z]i,j,0-7h4k2360[∂3u∂x2∂z+∂3u∂y2∂z]i,j,0-h4180[∂5u∂x4∂z+∂5u∂y4∂z]i,j,0+O(h6).Finally, to complete the explicit scheme for the boundary conditions, we present the previous equation in the form(34)δzui,j,0+h26(1+h2k230)[δzδx2ui,j,0+δzδy2ui,j,0]+h430[δzδx2δy2ui,j,0]=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0+7h4360[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0-h4180[∂5u∂x2∂y2∂z]i,j,0-h4k2180[∂3u∂x2∂z+∂3u∂y2∂z]i,j,0-h4180[∂5u∂x4∂z+∂5u∂y2∂z]i,j,0+O(h6).
Now by using (3) and replacing ui,j,k with Ui,j,k, we can write
(35)Ui,j,-1+h26(1+h2k230)[δx2Ui,j,-1+δy2Ui,j,-1]+h430δx2δy2Ui,j,-1=Ui,j,1+h26(1+h2k230)[δx2Ui,j,1+δy2Ui,j,1]+h430δx2δy2Ui,j,1-2h(1-h2k26+h4k4120)βi,j,0-h33(1-h2k220)×[∂f∂z]i,j,0-h560[∂3f∂z3]i,j,0+7h5180[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0+h590[∂4β∂x2∂y2+k2(∂2β∂x2+∂2β∂y2)+∂4β∂x4+∂4β∂y4]i,j,0.This formula allows us to eliminate the termUi,j,-1 in (28). The explicit implementation of the Neumann boundary conditions on the other boundary pieces can be conducted in a similar way.Sommerfeld-Like Boundary Conditions. In this section we consider the implementation of the Sommerfeld-like boundary conditions (2) which are often used in scattering problems
(36)∂u∂z+iku|z=0=0.
The difficulty in using similar methods to those described in the previous section can be seen in (35). The last term in this equation requires calculation of the fourth derivative of the right-hand side of (29) and if it depends on an unknown variable u, the use of a compact sixth order approximation scheme becomes problematic. To avoid such a difficulty, we introduce a new variable v=eikzu. Then the Helmholtz equation becomes
(37)∂2v∂x2+∂2v∂y2+∂2v∂z2-2ik∂v∂z=f¯,
where f¯=eikzf. Now (36) can be rewritten in the form
(38)∂v∂z|z=0=0.
The sixth order approximation of (38) becomes
(39)δzvi,j,0=∂v∂z|i,j,0+h26∂3v∂z3|i,j,0+h4120∂5v∂z5|i,j,0+O(h6).
Assuming sufficient smoothness of v and f¯ and using (37) and (38), we can express the third and fifth derivatives of v in the form
(40)∂3v∂z3|z=0=[∂f¯∂z+2ik∂2v∂z2]z=0,(41)∂5v∂z5|z=0=-4ik[2k2∂2v∂z2+∂4v∂z2∂x2+∂4v∂z2∂y2]z=0-[4k2∂f¯∂z-2ik∂2f¯∂z2+∂3f¯∂z∂x2hhhhh+∂3f¯∂z∂y2-∂3f¯∂z3]z=0.To preserve the sixth order compact approximation in (39) we need to use the fourth order approximation for ∂2v/∂z2 in (40). The second order approximation is sufficient for the derivatives of v in (41). First, let us consider the fourth finite-difference compact approximation for ∂2v/∂z2. It can be presented in the form
(42)∂2v∂z2|i,j,0=δz2vi,j,0-h212∂4v∂z4|i,j,0+O(h4)=δz2vi,j,0+h212[4k2∂2v∂z2+∂4v∂z2∂x2+∂4v∂z2∂y2δz2vi,j,0+h212hhh-2ik∂f¯∂z-∂2f¯∂z2]i,j,0+O(h4)=δz2vi,j,0+h212[(2ik∂f¯∂z+∂2f¯∂z2)i,j,04k2δz2vi,j,0+δz2δx2vi,j,0+δz2δy2vi,j,0δz2vi,j,0+h212hhh-(2ik∂f¯∂z+∂2f¯∂z2)i,j,0]+O(h4).Now we can use (38), (42) and the second order approximation of (41) to express (39) in the form
(43)δzvi,j,0-ikh23(1+2k2h215)δz2vi,j,0+ikh4180(δz2δx2vi,j,0+δz2δy2vi,j,0)=F¯i,j,0,F¯i,j,0=(h26+k2h445)∂f¯∂z|i,j,0-ikh490∂2f¯∂z2|i,j,0+h4120[∂3f¯∂z3-∂3f¯∂z∂x2-∂3f¯∂z∂y2]i,j,0.This equation gives an implicit compact sixth order approximation for (38). Next, the explicit implementation of the boundary conditions (38) similar to (35) is developed by using
(44)0=h2μ1[∂3v∂z∂x2+∂3v∂z∂y2]i,j,0=μ1h2[(∂5v∂z3∂x2+∂5v∂3z∂y2)i,j,0δzδx2vi,j,0+δzδy2vi,j,0μ1h2hhh-h26(∂5v∂z3∂x2+∂5v∂3z∂y2)i,j,0]+O(h6)=μ1h2[δzδx2vi,j,0+δzδy2vi,j,0]-ikμ1h43[∂4v∂z2∂x2+∂4v∂2z∂y2]i,j,0+μ1h46[2μ2∂5v∂x2∂y2∂z-∂3f¯∂z∂x2-∂3f¯∂z∂y2]i,j,0+O(h6).Here,μ1 and μ2 are parameters. We also use (∂5v/∂z∂x2∂y2)|z=0=(∂5v/∂z∂x4)|z=0=0 and (∂5v/∂z∂y4)|z=0=0. We leave the first of these derivatives in the right-hand side but drop the latter two derivatives. Now, by using the second order finite-difference approximation in (44) and adding it to the (43), we derive the formula
(45)μ3[δzvi,j,0-ikh23(1+2k2h215)δz2vi,j,0]+μ3ikh4180[δz2δx2vi,j,0+δz2δx2vi,j,0]+μ1μ3[ikh43(δz2δx2vi,j,0+δz2δx2vi,j,0)δzδx2vi,j,0+δzδx2vi,j,0μ1μ3hhh-ikh43(δz2δx2vi,j,0+δz2δx2vi,j,0)]+μ1μ2μ3δzδx2δzvi,j,0=μ3F¯+μ1μ3h46[∂3f¯∂z∂x2+∂3f¯∂z∂y2]i,j,0.Next we choose parametersμ1, μ2, and μ3 to match the coefficients in (28). This choice is determined by the conditions
(46)μ3=11+(2ikh/3)(1+(2k2h2/15)),μ1=(ikh/90)+(1/6μ3)(1+(k2h2/30))1+(2ikh/3),μ2=110μ1μ3.Using these parameters and replacingvi,j,l with Ui,j,leikzl, we write
(47)Ui,j,-1+h26(1+h2k230)(δx2+δy2)Ui,j,-1+h430δx2δy2Ui,j,-1=e2ikhμ3[h2(μ1+2ikh3(160-μ1))(δx2+δy2)Ui,j,1e2ikhμ3hhh+h430μ3δx2δy2Ui,j,1e2ikhμ3hhh+(1-2ikh3-4ih3k345)Ui,j,1]+eikhμ3[4ikh3(1+2h2k215)Ui,j,0+eikhμ3h-4ikh33(160-μ1)(δx2+δy2)Ui,j,0]-eikhμ3[2hF¯i,j,0+μ1h53(∂3f¯∂z∂x2+∂3f¯∂z∂y2)i,j,0].This equation provides an explicit compact sixth order approximation for the boundary condition (36). At the upper boundary z=a, the Sommerfeld-like boundary condition can be written as (∂u/∂z)-iku|z=a=0 and the auxiliary substitution becomes v=e-ikzu. The rest of the derivation is very similar to the derivation for the case of the lower boundary. Implementations of the sixth order compact approximation of the boundary conditions in the other directions are also similar to the calculations already presented.


### 3.3. Fast Fourier Preconditioner

The main idea of this paper is to utilize the existing lower-order approximation direct solvers as preconditioners in the implementation of a higher resolution scheme. One of the most popular methods for the approximate solution of the Helmholtz equation is the second order central difference scheme. We consider the cases with Dirichlet or Neumann boundary conditions on the sides of the computational domainΩ; that is, u(0,y,z)=u(a,y,z)=u(x,0,z)=u(x,a,z)=0 or (∂u/∂x)|x=0=(∂u/∂x)|x=a=(∂u/∂y)|y=0=(∂u/∂y)|y=a=0. At the top z=a and the bottom z=0 of the computational domain, we assign one of the three boundary conditions under consideration: Dirichlet, Neumann, or Sommerfeld-like boundary conditions (2). In the case of Neumann and Sommerfeld-like boundary conditions, we use the second order central finite-difference approximation of the first derivative on all boundaries. To express the preconditioner in general form we introduce the N×N matrices
(48)ΛNβ=[01+β⋯0101⋮⋱⋱⋱⋮1010⋯1+β0]
and DNα=diag⁡(α,0,…,0,α), where α, β, and N are parameters depending on boundary conditions. Next let Aβ,N0=ΛNβ-2I and Aα,β,N=ΛNβ-(2-k2h2)I+DNα. Using Kronecker products, the preconditioning matrix can now be written in the form
(49)Ap=Aβ,N0⊗I⊗I+I⊗Aβ,N0⊗I+I⊗I⊗Aαz,βz,Nz,
where β=0 in the case of Dirichlet boundary conditions at the sides of the computational domain and β=1 in the case of Neumann boundary conditions on the same boundaries. The number of grid points in the z-direction Nz is chosen in such a way that the grid step size is the same in all three directions. This condition is imposed so that the lower-order preconditioner and the sixth order approximation compact scheme developed in the previous section use the same grid points. Depending on the boundary conditions at the top and the bottom of the computational domain, the parameters αz and βz are defined in the following way: in the case of Dirichlet boundary conditions αz=0 and βz=0; in the case of Neumann boundary condition αz=0 and βz=1; and in the case of the Sommerfeld-like boundary conditions αz=2ikh and βz=1. At each step of the iterative process, the solution to the second equation in (12) with the matrix (49) can be obtained by using FFT-type algorithms in O(N2Nzlog⁡2N) operations.


### 3.4.  Convergence Analysis of 3D Algorithm

In this section, the convergence of the proposed algorithms is considered in the case of the 3D Helmholtz equation with the Dirichlet boundary conditions (1) and (2) on the rectangular computational domain Ω={0≤x,y,z≤1} with uniform grid size h=1/(N+1). To insure the uniqueness of the solution to the original boundary value problem, we assume that the coefficient k2 is in the resolvent set of the Laplace operator defined on an appropriate function space. We also impose a common restriction on the number of grid points per wavelength; that is, kh<2π/10. The convergence analysis of the proposed algorithms follows the same lines as in Theorem 5. The preconditioning matrix is given by (49) with β=αz=βz=0 and Nz=N. To express the right-hand side of the sixth order approximation scheme given in (28), we will use the notation introduced in (49). In addition, we define the N×N matrix Aβ,N1=ΛNβ-4I. Then the resulting N3×N3 matrix A can be written as
(50)A=A0,N0⊗I⊗I+I⊗A0,N0⊗I+I⊗I⊗A0,N0+16[1+h2k230]×(ΛN0⊗A0,N1⊗I+I⊗A0,N1⊗A0,N1-4I⊗I⊗I)+130A0,N1⊗A0,N1⊗A0,N1+[h2k2-h4k412+h6k6360]I⊗I⊗I.Both matricesA and Ap have the same set of orthonormal eigenvectors Vi,j,l (see, e.g., [17]):
(51)Vi,j,lm,n,s=(2h)3/2sin(imπh)sin(jnπh)sin(lsπh),Vi,j,lm,n,s=(2h)3/2sin(imπh)1≤i,j,l,m,n,s≤N.The matrix of eigenvectors is unitary so thel2-norm of V and V-1 is one. To present the eigenvalues of the matrices, we introduce the notation sr=sin2(rπh/2),r=i,j,l. Then the eigenvalues of the matrices Ap and A in the 3D case are given by
(52)λi,j,lp=-4[si2+sj2+sl2]+h2k2,λi,j,l=-4[si2+sj2+sl2]+h2k2+83[1+h2k230][si2sj2+si2sl2+sl2sj2]-3215si2sj2sl2-h4k412+h6k6360,i,j,l=1,…,N,respectively.Unfortunately, the eigenvalues of these matrices do not satisfy the hypotheses of Theorem2, but we still can apply Corollary 3 to prove the convergence of the GMRES method and the SKS algorithm by showing that λi,j,l/λi,j,lp=1+di,j,l, where |di,j,l|<1, for all i,j,l=1,…,N in some simple cases.Let us assume thatN≥1; that is, h≤1/2 and k2<25. Then sr2=sin2(πrh/2)≥sin2(πh/2)≥2h2, for r=i,j,l. It is also easy to see that (8/3)si2sj2≥(32/15)si2sj2sl2 and (8/3)(si2sl2+sj2sl2)≥64h4/3≥k4h4/12 for the given values of k. Now the eigenvalues of the preconditioned matrix AAp-1 can be written in the form
(53)λi,j,lλi,j,lp=1+[83(1+h2k230)[si2sj2+si2sl2+sl2sj2]hhhhhh-3215si2sj2sl2-h4k412+h6k6360]/λi,j,lp,0000000000000000000000i,j,l=1,…,N.This gives us a presentation of the eigenvalues of the preconditioned matrix in the formλi,j,l/λi,j,lp=1+di,j,l, where |di,j,l|<3/4, for all i,j,l=1,…,N. Now we can apply Corollary 3 and obtain an estimate for the rate of convergence of the GMRES and SKS algorithms:
(54)∥rn∥2≤(34)n∥r0∥2,
where rn=F-AU(n).Similar to the 1D case, it is possible to show that this estimate holds fork>6 and a sequence of grids with the grid sizes h1>h2>h3⋯ if
(55)min⁡i,j,l=1,…,N|4(si2+sj2+sl2)h2-k2|≥δ0>0
for some δ0=const and sufficiently small h1. As in the case of the 1D problem if follows from the fact that the resolvent set of a bounded linear operator is open. The proof is somewhat tedious but elementary, so we skip it. In Table 1 we present lower and upper bounds m and M for λi,j,l/λi,j,lp-1 for different values of k and for different grid step sizes. Also, the last two rows of Table 1 possible choices for h1 and δ0 for different k.Table 1
3D Dirichlet problem eigenvalue bounds (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
-
0.015
∣
0.49
-
1.22
∣
0.49
-
1.02
∣
0.88
-
34.7
∣
180
-
10.2
∣
55.0
h
=
1
/
128
-
0.0038
∣
0.49
-
0.158
∣
0.49
-
0.140
∣
0.64
-
2.81
∣
0.90
-
2.21
∣
1.74
h
=
1
/
256
-
0.00094
∣
0.49
-
0.035
∣
0.49
-
0.031
∣
0.49
-
0.22
∣
0.69
-
11.7
∣
0.56
h
=
1
/
512
-
0.00023
∣
0.49
-
0.009
∣
0.49
-
0.008
∣
0.49
-
0.048
∣
0.49
-
0.298
∣
0.49
δ
0
8.00
4.00
1.80
1.12
2.9
h
1
1/64
1/128
1/256
1/256
1/512This result is weaker than in the 1D case; that is, it says that convergence is independent of the grid step, but it does not improve with the decrease ofh. In this case, Ap is said to be “an optimal preconditioner” (see, e.g., [21, page 196]). This situation is typical for multigrid-type algorithms but it is not what we would expect based on 1D example. In the following section we present the results of numerical experiments in the case of several test problems which confirm that actual convergence of the presented methods accelerate with the decrease of the grid size; that is, the convergence in the numerical experiments significantly exceeds the estimate (54).


## 3.1. A Sixth Order Approximation Compact Scheme

In this section we present a three-dimensional compact sixth order approximation finite-difference scheme that was first introduced in [6], and we develop a sixth order compact explicit approximation of the Neumann and Sommerfeld-like boundary conditions (2). In this discussion we consider a uniform grid Ωh={(xi,yj,zk)∣xi=ih,yj=jh,z=kh;i,j,k=0,…,N-1,h=a/(N-1)}.First, we consider a finite-difference approximation of (1) in the form
(26)δx2ui,j,k+δy2ui,j,k+δz2ui,j,k+k2ui,j,k+Ti,j,k=fi,j,k,
where ui,j,k=u(xi,yj,zk) and fi,j,k=f(xi,yj,zk) and
(27)Ti,j,k=-[h212(∂4u∂x4+∂4u∂y4+∂4u∂z4)hhh+h4360(∂6u∂x6+∂6u∂y6+∂6u∂z6)]i,j,k+O(h6).Using the appropriate derivatives of (1) we can write the sixth order compact approximation of the Helmholtz equation in the form
(28)(δx2+δy2+δz2)Ui,j,k+h26(1+k2h230)×(δx2δy2+δx2δz2+δy2δz2)Ui,j,k+h430δx2δy2δz2Ui,j,k+k2(1-k2h212+k4h4360)Ui,j,k=(1-k2h212+k4h4360)fi,j,k+h212(1-k2h230)∇2fi,j,k+h4360∇4fi,j,k+h490(∂4f∂x2∂y2+∂4f∂x2∂z2+∂4f∂y2∂z2)i,j,k.
For the two-dimensional case this scheme was proposed in [7]. In the three-dimensional case with Dirichlet boundary conditions, it was developed in [6]. We will consider the iterative implementation of this scheme based on preconditioned Krylov subspace-type algorithms.


## 3.2. Boundary Conditions

The implementation of the Dirichlet boundary conditions (2) is straightforward but the explicit compact approximation of Neumann and Sommerfeld-like boundary conditions requires careful consideration. The sixth order compact approximation of the Neumann boundary condition in the two-dimensional case was considered in [7]. Here we extend this approach to the three-dimensional problem.Neumann Boundary Conditions. In this subsection the Neumann boundary conditions are considered in the form
(29)∂u∂z|z=0=β(x,y).We restrict our consideration to only one side of the computational domainΩ. By using Taylor series, we can derive the sixth order approximation formula
(30)δzui,j,0=∂u∂z|i,j,0+h26∂3u∂z3|i,j,0+h4120∂5u∂z5|i,j,0+O(h6).
To express the third and fifth derivatives in (30) we can differentiate the original Helmholtz equation (1), assuming sufficient smoothness of the solution and the right-hand side. After substitution of these derivatives into (30), the resulting expression is
(31)δzui,j,0=[∂u∂z]i,j,0+h26×[∂f∂z-∂3u∂x2∂z-∂3u∂y2∂z-k2∂u∂z]i,j,0+h4120[∂3f∂z3-∂5u∂z3∂x2-∂5u∂z3∂y2hhhhhhh-k2(∂f∂z-∂3u∂x2∂z-∂3u∂y2∂z-k2∂u∂z)]i,j,0+O(h6).Next, we approximate the third order mixed derivatives in the second bracket of (31) by a fourth order approximation formula and in the second line using a second order approximation scheme, which yields
(32)δzui,j,0=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0-h26[δzδx2u+δzδy2u-h26(∂5u∂x2∂z3+∂5u∂y2∂z3)hhhhhh-h212(∂5u∂x4∂z+∂5u∂y4∂z)]i,j,0+h4k2120[δzδx2ui,j,0+δzδy2ui,j,0]-h4120[∂5u∂x2∂z3+∂5u∂y2∂z3]i,j,0+O(h6).
Now we can simplify this and write
(33)δzui,j,0+h26(1-h2k220)[δzδx2ui,j,0+δzδy2ui,j,0]=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0+7h4360[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0-7h4180[∂5u∂x2∂y2∂z]i,j,0-7h4k2360[∂3u∂x2∂z+∂3u∂y2∂z]i,j,0-h4180[∂5u∂x4∂z+∂5u∂y4∂z]i,j,0+O(h6).Finally, to complete the explicit scheme for the boundary conditions, we present the previous equation in the form(34)δzui,j,0+h26(1+h2k230)[δzδx2ui,j,0+δzδy2ui,j,0]+h430[δzδx2δy2ui,j,0]=(1-h2k26+h4k4120)[∂u∂z]i,j,0+h26(1-h2k220)[∂f∂z]i,j,0+h4120[∂3f∂z3]i,j,0+7h4360[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0-h4180[∂5u∂x2∂y2∂z]i,j,0-h4k2180[∂3u∂x2∂z+∂3u∂y2∂z]i,j,0-h4180[∂5u∂x4∂z+∂5u∂y2∂z]i,j,0+O(h6).
Now by using (3) and replacing ui,j,k with Ui,j,k, we can write
(35)Ui,j,-1+h26(1+h2k230)[δx2Ui,j,-1+δy2Ui,j,-1]+h430δx2δy2Ui,j,-1=Ui,j,1+h26(1+h2k230)[δx2Ui,j,1+δy2Ui,j,1]+h430δx2δy2Ui,j,1-2h(1-h2k26+h4k4120)βi,j,0-h33(1-h2k220)×[∂f∂z]i,j,0-h560[∂3f∂z3]i,j,0+7h5180[∂3f∂z∂2x+∂3f∂z∂2y]i,j,0+h590[∂4β∂x2∂y2+k2(∂2β∂x2+∂2β∂y2)+∂4β∂x4+∂4β∂y4]i,j,0.This formula allows us to eliminate the termUi,j,-1 in (28). The explicit implementation of the Neumann boundary conditions on the other boundary pieces can be conducted in a similar way.Sommerfeld-Like Boundary Conditions. In this section we consider the implementation of the Sommerfeld-like boundary conditions (2) which are often used in scattering problems
(36)∂u∂z+iku|z=0=0.
The difficulty in using similar methods to those described in the previous section can be seen in (35). The last term in this equation requires calculation of the fourth derivative of the right-hand side of (29) and if it depends on an unknown variable u, the use of a compact sixth order approximation scheme becomes problematic. To avoid such a difficulty, we introduce a new variable v=eikzu. Then the Helmholtz equation becomes
(37)∂2v∂x2+∂2v∂y2+∂2v∂z2-2ik∂v∂z=f¯,
where f¯=eikzf. Now (36) can be rewritten in the form
(38)∂v∂z|z=0=0.
The sixth order approximation of (38) becomes
(39)δzvi,j,0=∂v∂z|i,j,0+h26∂3v∂z3|i,j,0+h4120∂5v∂z5|i,j,0+O(h6).
Assuming sufficient smoothness of v and f¯ and using (37) and (38), we can express the third and fifth derivatives of v in the form
(40)∂3v∂z3|z=0=[∂f¯∂z+2ik∂2v∂z2]z=0,(41)∂5v∂z5|z=0=-4ik[2k2∂2v∂z2+∂4v∂z2∂x2+∂4v∂z2∂y2]z=0-[4k2∂f¯∂z-2ik∂2f¯∂z2+∂3f¯∂z∂x2hhhhh+∂3f¯∂z∂y2-∂3f¯∂z3]z=0.To preserve the sixth order compact approximation in (39) we need to use the fourth order approximation for ∂2v/∂z2 in (40). The second order approximation is sufficient for the derivatives of v in (41). First, let us consider the fourth finite-difference compact approximation for ∂2v/∂z2. It can be presented in the form
(42)∂2v∂z2|i,j,0=δz2vi,j,0-h212∂4v∂z4|i,j,0+O(h4)=δz2vi,j,0+h212[4k2∂2v∂z2+∂4v∂z2∂x2+∂4v∂z2∂y2δz2vi,j,0+h212hhh-2ik∂f¯∂z-∂2f¯∂z2]i,j,0+O(h4)=δz2vi,j,0+h212[(2ik∂f¯∂z+∂2f¯∂z2)i,j,04k2δz2vi,j,0+δz2δx2vi,j,0+δz2δy2vi,j,0δz2vi,j,0+h212hhh-(2ik∂f¯∂z+∂2f¯∂z2)i,j,0]+O(h4).Now we can use (38), (42) and the second order approximation of (41) to express (39) in the form
(43)δzvi,j,0-ikh23(1+2k2h215)δz2vi,j,0+ikh4180(δz2δx2vi,j,0+δz2δy2vi,j,0)=F¯i,j,0,F¯i,j,0=(h26+k2h445)∂f¯∂z|i,j,0-ikh490∂2f¯∂z2|i,j,0+h4120[∂3f¯∂z3-∂3f¯∂z∂x2-∂3f¯∂z∂y2]i,j,0.This equation gives an implicit compact sixth order approximation for (38). Next, the explicit implementation of the boundary conditions (38) similar to (35) is developed by using
(44)0=h2μ1[∂3v∂z∂x2+∂3v∂z∂y2]i,j,0=μ1h2[(∂5v∂z3∂x2+∂5v∂3z∂y2)i,j,0δzδx2vi,j,0+δzδy2vi,j,0μ1h2hhh-h26(∂5v∂z3∂x2+∂5v∂3z∂y2)i,j,0]+O(h6)=μ1h2[δzδx2vi,j,0+δzδy2vi,j,0]-ikμ1h43[∂4v∂z2∂x2+∂4v∂2z∂y2]i,j,0+μ1h46[2μ2∂5v∂x2∂y2∂z-∂3f¯∂z∂x2-∂3f¯∂z∂y2]i,j,0+O(h6).Here,μ1 and μ2 are parameters. We also use (∂5v/∂z∂x2∂y2)|z=0=(∂5v/∂z∂x4)|z=0=0 and (∂5v/∂z∂y4)|z=0=0. We leave the first of these derivatives in the right-hand side but drop the latter two derivatives. Now, by using the second order finite-difference approximation in (44) and adding it to the (43), we derive the formula
(45)μ3[δzvi,j,0-ikh23(1+2k2h215)δz2vi,j,0]+μ3ikh4180[δz2δx2vi,j,0+δz2δx2vi,j,0]+μ1μ3[ikh43(δz2δx2vi,j,0+δz2δx2vi,j,0)δzδx2vi,j,0+δzδx2vi,j,0μ1μ3hhh-ikh43(δz2δx2vi,j,0+δz2δx2vi,j,0)]+μ1μ2μ3δzδx2δzvi,j,0=μ3F¯+μ1μ3h46[∂3f¯∂z∂x2+∂3f¯∂z∂y2]i,j,0.Next we choose parametersμ1, μ2, and μ3 to match the coefficients in (28). This choice is determined by the conditions
(46)μ3=11+(2ikh/3)(1+(2k2h2/15)),μ1=(ikh/90)+(1/6μ3)(1+(k2h2/30))1+(2ikh/3),μ2=110μ1μ3.Using these parameters and replacingvi,j,l with Ui,j,leikzl, we write
(47)Ui,j,-1+h26(1+h2k230)(δx2+δy2)Ui,j,-1+h430δx2δy2Ui,j,-1=e2ikhμ3[h2(μ1+2ikh3(160-μ1))(δx2+δy2)Ui,j,1e2ikhμ3hhh+h430μ3δx2δy2Ui,j,1e2ikhμ3hhh+(1-2ikh3-4ih3k345)Ui,j,1]+eikhμ3[4ikh3(1+2h2k215)Ui,j,0+eikhμ3h-4ikh33(160-μ1)(δx2+δy2)Ui,j,0]-eikhμ3[2hF¯i,j,0+μ1h53(∂3f¯∂z∂x2+∂3f¯∂z∂y2)i,j,0].This equation provides an explicit compact sixth order approximation for the boundary condition (36). At the upper boundary z=a, the Sommerfeld-like boundary condition can be written as (∂u/∂z)-iku|z=a=0 and the auxiliary substitution becomes v=e-ikzu. The rest of the derivation is very similar to the derivation for the case of the lower boundary. Implementations of the sixth order compact approximation of the boundary conditions in the other directions are also similar to the calculations already presented.


## 3.3. Fast Fourier Preconditioner

The main idea of this paper is to utilize the existing lower-order approximation direct solvers as preconditioners in the implementation of a higher resolution scheme. One of the most popular methods for the approximate solution of the Helmholtz equation is the second order central difference scheme. We consider the cases with Dirichlet or Neumann boundary conditions on the sides of the computational domainΩ; that is, u(0,y,z)=u(a,y,z)=u(x,0,z)=u(x,a,z)=0 or (∂u/∂x)|x=0=(∂u/∂x)|x=a=(∂u/∂y)|y=0=(∂u/∂y)|y=a=0. At the top z=a and the bottom z=0 of the computational domain, we assign one of the three boundary conditions under consideration: Dirichlet, Neumann, or Sommerfeld-like boundary conditions (2). In the case of Neumann and Sommerfeld-like boundary conditions, we use the second order central finite-difference approximation of the first derivative on all boundaries. To express the preconditioner in general form we introduce the N×N matrices
(48)ΛNβ=[01+β⋯0101⋮⋱⋱⋱⋮1010⋯1+β0]
and DNα=diag⁡(α,0,…,0,α), where α, β, and N are parameters depending on boundary conditions. Next let Aβ,N0=ΛNβ-2I and Aα,β,N=ΛNβ-(2-k2h2)I+DNα. Using Kronecker products, the preconditioning matrix can now be written in the form
(49)Ap=Aβ,N0⊗I⊗I+I⊗Aβ,N0⊗I+I⊗I⊗Aαz,βz,Nz,
where β=0 in the case of Dirichlet boundary conditions at the sides of the computational domain and β=1 in the case of Neumann boundary conditions on the same boundaries. The number of grid points in the z-direction Nz is chosen in such a way that the grid step size is the same in all three directions. This condition is imposed so that the lower-order preconditioner and the sixth order approximation compact scheme developed in the previous section use the same grid points. Depending on the boundary conditions at the top and the bottom of the computational domain, the parameters αz and βz are defined in the following way: in the case of Dirichlet boundary conditions αz=0 and βz=0; in the case of Neumann boundary condition αz=0 and βz=1; and in the case of the Sommerfeld-like boundary conditions αz=2ikh and βz=1. At each step of the iterative process, the solution to the second equation in (12) with the matrix (49) can be obtained by using FFT-type algorithms in O(N2Nzlog⁡2N) operations.


## 3.4.  Convergence Analysis of 3D Algorithm

In this section, the convergence of the proposed algorithms is considered in the case of the 3D Helmholtz equation with the Dirichlet boundary conditions (1) and (2) on the rectangular computational domain Ω={0≤x,y,z≤1} with uniform grid size h=1/(N+1). To insure the uniqueness of the solution to the original boundary value problem, we assume that the coefficient k2 is in the resolvent set of the Laplace operator defined on an appropriate function space. We also impose a common restriction on the number of grid points per wavelength; that is, kh<2π/10. The convergence analysis of the proposed algorithms follows the same lines as in Theorem 5. The preconditioning matrix is given by (49) with β=αz=βz=0 and Nz=N. To express the right-hand side of the sixth order approximation scheme given in (28), we will use the notation introduced in (49). In addition, we define the N×N matrix Aβ,N1=ΛNβ-4I. Then the resulting N3×N3 matrix A can be written as
(50)A=A0,N0⊗I⊗I+I⊗A0,N0⊗I+I⊗I⊗A0,N0+16[1+h2k230]×(ΛN0⊗A0,N1⊗I+I⊗A0,N1⊗A0,N1-4I⊗I⊗I)+130A0,N1⊗A0,N1⊗A0,N1+[h2k2-h4k412+h6k6360]I⊗I⊗I.Both matricesA and Ap have the same set of orthonormal eigenvectors Vi,j,l (see, e.g., [17]):
(51)Vi,j,lm,n,s=(2h)3/2sin(imπh)sin(jnπh)sin(lsπh),Vi,j,lm,n,s=(2h)3/2sin(imπh)1≤i,j,l,m,n,s≤N.The matrix of eigenvectors is unitary so thel2-norm of V and V-1 is one. To present the eigenvalues of the matrices, we introduce the notation sr=sin2(rπh/2),r=i,j,l. Then the eigenvalues of the matrices Ap and A in the 3D case are given by
(52)λi,j,lp=-4[si2+sj2+sl2]+h2k2,λi,j,l=-4[si2+sj2+sl2]+h2k2+83[1+h2k230][si2sj2+si2sl2+sl2sj2]-3215si2sj2sl2-h4k412+h6k6360,i,j,l=1,…,N,respectively.Unfortunately, the eigenvalues of these matrices do not satisfy the hypotheses of Theorem2, but we still can apply Corollary 3 to prove the convergence of the GMRES method and the SKS algorithm by showing that λi,j,l/λi,j,lp=1+di,j,l, where |di,j,l|<1, for all i,j,l=1,…,N in some simple cases.Let us assume thatN≥1; that is, h≤1/2 and k2<25. Then sr2=sin2(πrh/2)≥sin2(πh/2)≥2h2, for r=i,j,l. It is also easy to see that (8/3)si2sj2≥(32/15)si2sj2sl2 and (8/3)(si2sl2+sj2sl2)≥64h4/3≥k4h4/12 for the given values of k. Now the eigenvalues of the preconditioned matrix AAp-1 can be written in the form
(53)λi,j,lλi,j,lp=1+[83(1+h2k230)[si2sj2+si2sl2+sl2sj2]hhhhhh-3215si2sj2sl2-h4k412+h6k6360]/λi,j,lp,0000000000000000000000i,j,l=1,…,N.This gives us a presentation of the eigenvalues of the preconditioned matrix in the formλi,j,l/λi,j,lp=1+di,j,l, where |di,j,l|<3/4, for all i,j,l=1,…,N. Now we can apply Corollary 3 and obtain an estimate for the rate of convergence of the GMRES and SKS algorithms:
(54)∥rn∥2≤(34)n∥r0∥2,
where rn=F-AU(n).Similar to the 1D case, it is possible to show that this estimate holds fork>6 and a sequence of grids with the grid sizes h1>h2>h3⋯ if
(55)min⁡i,j,l=1,…,N|4(si2+sj2+sl2)h2-k2|≥δ0>0
for some δ0=const and sufficiently small h1. As in the case of the 1D problem if follows from the fact that the resolvent set of a bounded linear operator is open. The proof is somewhat tedious but elementary, so we skip it. In Table 1 we present lower and upper bounds m and M for λi,j,l/λi,j,lp-1 for different values of k and for different grid step sizes. Also, the last two rows of Table 1 possible choices for h1 and δ0 for different k.Table 1
3D Dirichlet problem eigenvalue bounds (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
-
0.015
∣
0.49
-
1.22
∣
0.49
-
1.02
∣
0.88
-
34.7
∣
180
-
10.2
∣
55.0
h
=
1
/
128
-
0.0038
∣
0.49
-
0.158
∣
0.49
-
0.140
∣
0.64
-
2.81
∣
0.90
-
2.21
∣
1.74
h
=
1
/
256
-
0.00094
∣
0.49
-
0.035
∣
0.49
-
0.031
∣
0.49
-
0.22
∣
0.69
-
11.7
∣
0.56
h
=
1
/
512
-
0.00023
∣
0.49
-
0.009
∣
0.49
-
0.008
∣
0.49
-
0.048
∣
0.49
-
0.298
∣
0.49
δ
0
8.00
4.00
1.80
1.12
2.9
h
1
1/64
1/128
1/256
1/256
1/512This result is weaker than in the 1D case; that is, it says that convergence is independent of the grid step, but it does not improve with the decrease ofh. In this case, Ap is said to be “an optimal preconditioner” (see, e.g., [21, page 196]). This situation is typical for multigrid-type algorithms but it is not what we would expect based on 1D example. In the following section we present the results of numerical experiments in the case of several test problems which confirm that actual convergence of the presented methods accelerate with the decrease of the grid size; that is, the convergence in the numerical experiments significantly exceeds the estimate (54).


## 4. Numerical Results

In this section we present the results of numerical experiments which demonstrate the quality of the proposed numerical methods. These algorithms were implemented in Matlab 7.11.0 on an iMac with an Intel Core i7, 2.93 GHz processor. We also use the standard programing implementation of the GMRES method in Matlab (gmres function).


### 4.1. 1D Test Problem

In the first series of numerical experiments, we consider the convergence of the compact sixth order approximation algorithm on a sequence of grids. In these tests we focus on the numerical solution of the 1D Dirichlet problem for the Helmholtz equation with zero boundary conditions on the interval0≤x≤1 and k=20. The computational grid in our experiments is defined by xi=ih,i=0,…,N+1, where h=1/(N+1) and N is the number of grid points. We consider the function u(xi)=xi(1-xi)cos⁡(kπxi),i=1,…,N as the exact solution of the original boundary value problem with the right-hand side f(xi)=-(2+k2(π2-xi(1-xi))cos⁡(kπxi) + 2kπ(2xi-1)sin(kπxi),i=1,…,N. The first column of Table 2 displays the step size h of the grids in our experiments. Columns 2–4 of the table present the number of iterations until the convergence required by the preconditioned GMRES method, the SKS method (Algorithm 1), and the Chebyshev acceleration (CA) algorithm [24]. In all our experiments we use restarted GMRES method with the restart value 20. We use the stoppage criterion ∥F-AU(n)∥2/∥F∥2≤10-10, where U(n) is the iterative solution on nth iteration. Column 5 reports the approximation error of the numerical solution Err6=max⁡1≤j≤N|Uj(n)-u(xj)|. In Column 6, the second order approximation error Err2 of the numerical solution obtained by using the preconditioner as a solver is presented for comparison. The number of grid points is chosen to satisfy the requirement h<2π/10k. The values of the minimum and maximum eigenvalues of the preconditioned matrix used by the Chebyshev acceleration algorithm are calculated exactly by (22).Table 2
1D Dirichlet problem.
Grid size (h)
GMRES
SKS
CA
Err
6
Err
2
1
/
32
6
8
6
3.88
*
1
0
-
03
1.16
*
1
0
-
01
1
/
64
4
4
5
5.08
*
1
0
-
05
2.4
*
1
0
-
02
1
/
128
3
3
4
7.61
*
1
0
-
7
5.7
*
1
0
-
03
1
/
258
3
3
3
1.14
*
1
0
-
8
1.4
*
1
0
-
03These experiments demonstrate the sixth order convergence of the numerical solution to the exact solution of the original boundary value problem. Theorem5 suggests that the residual on nth iteration in the presented iterative algorithms decreases with the increase in the number of grid points. The results of the numerical experiments confirm this conclusion. We have already shown that the order of the preconditioned matrix in this case is two. The numerical value of the constant ψ calculated by using (25) is 1.99. The 1D numerical experiments confirm the effectiveness of the iterative approaches under consideration and allow us to expect that the same properties hold for 3D problems.


### 4.2. 3D Test Problems

In the next series of experiments we consider different boundary value problems for the 3D Helmholtz equations (1) and (2). The computational domain is given by Ω={(x,y,z)∣0≤x,y,z≤1} with a uniform grid in all three directions xi=ih,i=1,…,Nx,yj=jh,j=1,…,Ny and zk=kh,k=1,…,Nz. The numbers of grid points Nx, Ny, and Nz in x,y, and z-directions are chosen to provide uniform grid size.Dirichlet Problem. First, we consider the Helmholtz equation with boundary conditions u=0 and k=20. This test is similar to the 1D case with Nx=Ny=Nz=N and h=1/(N+1). As the solution of the boundary value problem we choose the function u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy),ϕ(z)=sin(kπz), and 0≤x,y,z≤1.In the first series of tests we consider the convergence of the compact high order approximation finite-difference scheme (28) on a sequence of grids as well as convergence properties of the GMRES, SKS, and CA iterative algorithms. The iterative process is stopped when the initial residual is reduced by a factor of 10-10. The values of the minimum and maximum eigenvalues of the preconditioned matrix used by the Chebyshev acceleration algorithm are calculated exactly using (53). The description of the columns of Table 3 is similar to the description of Table 2.Table 3
3D Dirichlet problem.
Grid size (h)
GMRES
SKS
CA
Err
6
Err
2
1
/
64
5
10
23
3.03
*
1
0
-
04
4.47
*
1
0
-
06
1
/
128
3
6
15
8.13
*
1
0
-
05
6.35
*
1
0
-
08
1
/
256
3
5
14
2.06
*
1
0
-
05
9.68
*
1
0
-
10The test results confirm the sixth order approximation of the compact finite-difference scheme. Though, in our theoretical analysis (53) we could not prove the desired kth order preconditioning property in the case of the 3D Dirichlet problem, the SKS algorithm still exhibits acceleration of the convergence with the decrease of grid size. We also mention that the total computer time required for convergence for the SKS and GMRES methods is approximately the same and the Chebyshev acceleration algorithm requires about four times more computer time until convergence than the other two iterative algorithms do. The SKS algorithm has also greater potential for efficient implementation on parallel computers than the GMRES method does.One potential application of the proposed methods is the electromagnetic scattering problems. In these problems, the typical value of the coefficientk is in the range of 10 and 50. For example, the propagation of an electromagnetic wave with 1 GHz frequency in a vacuum is described by the Helmholtz equation with k≈20.9m-1. In Table 4 we present the convergence of the algorithms presented here on the sequence of grids with different values k. In every cell of the table we present the number of iteration for three methods under consideration: GMRES, SKS, and CA. The stoppage criterion is the same as in the previous series of experiments. If any of these three methods fails to converge in 100 iterations, we indicate this by using symbol “>100" and in the case of divergence we use “div” symbol.Table 4
3D Dirichlet problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
4
∣
7
∣
14
5
∣
10
∣
23
6
∣
13
∣
25
12
∣
div
∣ >100
19
∣
div
∣ >100
h
=
1
/
128
3
∣
5
∣
14
3
∣
6
∣
15
4
∣
8
∣
17
4
∣
10
∣
69
5
∣
11
∣ >100
h
=
1
/
256
3
∣
4
∣
12
3
∣
5
∣
14
3
∣
6
∣
13
3
∣
6
∣
22
4
∣
7
∣
63
ψ
2.02
2.00
1.98
2.00
2.02From these numerical experiments we can see that the preconditioned GMRES and SKS methods demonstrate excellent convergence properties. We also observe that the convergence of the proposed algorithms improves with the increase in the number of the grid points. In the analysis of the algorithm’s convergence (54), it was shown that the number of iterations does not increase with the increase of the grid points. But the last row in Table 4 indicates that the parameter ψ (25) for all values of k is approximately 2, which indicates that similar to the 1D case the convergence of the 3D algorithm not only does not depend on the grid size but accelerates with the increase of the number of grid points. This also indicates that even in the case when the SKS method does not converge, one can expect convergence on a finer grid. In all numerical tests, the parameter ψ is calculated by using the first two iterations of the SKS method on the grids with the grid sizes h=1/128 and h=1/256.From the results presented, it is clear that slow convergence and difficulties arising in the choice of the parametersm^ and M^ in the Chebyshev acceleration algorithm (18) make this method a poor choice for the implementation of the developed high-order approximation scheme. These numerical tests illustrate the fact that the Chebyshev polynomial is not optimal on the discrete spectrum and this has a dramatic effect on the iterative methods based on this polynomial. So, in the next numerical experiments we focus on the convergence properties of the GMRES and SKS methods in our numerical framework.Dirichlet-Neumann Boundary Conditions. In the next series of experiments we replace the Dirichlet boundary conditions u=0 in (2) with u=ub at the top of the computational domain (z=1) and with the Neumann boundary conditions ∂u/∂z=0 at z=0. On the other boundaries of the rectangular computational domain Ω={0≤x,y,z≤1} we still use the Dirichlet boundary conditions u=0. As a test function that satisfies the boundary conditions we consider u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy),ϕ(z)=cos⁡(kπz), and 0≤x,y,z≤1. Also, ub=ϕ1(x)ϕ2(y)ϕ3(1),0≤x,y≤1. The number of grid points is chosen to provide the same grid step in all three directions; that is Nx=Ny=Nz-1. The stoppage criterion is the same as in the previous numerical experiments. Table 5 provides the results of the numerical tests on a sequence of grids with the coefficient k=20. As already mentioned we focus only on the two best algorithms: GMRES and SKS. The rest of Table 5 is similar to Table 3.Table 5
3D Dirichlet-Neumann problem.
Grid size (h)
GMRES
SKS
Err
6
Err
2
1
/
64
8
11
4.65
*
1
0
-
06
3.67
*
1
0
-
04
1
/
128
5
8
6.61
*
1
0
-
08
8.74
*
1
0
-
05
1
/
256
4
7
1.01
*
1
0
-
09
2.17
*
1
0
-
05The results of this series of numerical tests indicate the sixth order convergence of the resulting compact finite-difference scheme including explicit approximation of the Neumann boundary conditions. As expected, the number of iterations decreases with the increase of the number of grid points in both the GMRES and SKS approaches. In Table6 we present the number of iterations required for convergence for both algorithms on a sequence of grids and with different k coefficients in the Helmholtz equation. The content of Table 6 is similar to the content of Table 4, except for the number of iterative approaches under consideration. We present only the number of iterations for the GMRES and SKS methods.Table 6
3D Dirichlet-Neumann problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
6
∣
10
8
∣
11
18
∣
div
36
∣
div
>600∣div
h
=
1
/
128
4
∣
9
5
∣
8
8
∣
9
14
∣
div
28
∣
div
h
=
1
/
256
3
∣
8
4
∣
7
5
∣
div
9
∣
div
7
∣
div
h
=
1
/
512
2
∣
6
3
∣
6
3
∣
6
5
∣
6
4
∣
6
ψ
2.57
2.18
2.39
1.65
1.94It follows from this table that the GMRES approach is much more efficient on coarse grids but on finer grids, the processor time for the SKS method is the same or even smaller than for the GMRES algorithm. For example, on the grid withh=1/512 and k=50 the corresponding processor times for the GMRES and SKS methods are 1467 sec and 1052 sec. Also, as we mentioned before, the SKS method has much greater potential than the GMRES algorithm for efficient implementation on parallel computers. We also expect that even if the SKS method does not converge on coarser grids, by reducing the grid size one can achieve convergence of this algorithm. We can see this in the example of the convergence of the methods under consideration in the case of k≥30. The values of the parameter ψ vary significantly for different values of k but ψ>0 still indicates that the convergence of the iterative algorithms improves with the increase in the number of the grid points.Dirichlet-Sommerfeld-Like Boundary Conditions. In the last series of numerical experiments, we consider the boundary value problem for the Helmholtz equation (1) with a combination of Sommerfeld-like (radiation) boundary conditions (2) at z=0 and z=1, and Dirichlet boundary conditions at all other boundaries of the rectangular computational domain Ω={0≤x,y,z≤1}. The source function f in (1) selected that the true solution is u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy), ϕ3(z)=eikz+e-ik(z-1)-2, and 0≤x,y,z≤1. Note that the analytic solution satisfies the radiation boundary condition (2). The numbers of grid points in x,y,z-directions are chosen such that Nx=Ny=Nz-2. First we consider the convergence of our sixth order approximation scheme on the sequence of grids with k=20. The main goal of this series of numerical experiments is to investigate the convergence properties of the new explicit compact sixth order approximation of the Sommerfeld-like boundary conditions (47) proposed in this paper. As in previous test runs, the iterative process is stopped when the l2-norm of the initial residual is reduced by a factor of 10-10. The results are presented in Table 7. The description of data presented in this table is the same as in the case of Table 5.Table 7
3D Dirichlet-Sommerfeld problem.
Grid size (h)
GMRES
SKS
Err
6
Err
2
1
/
64
6
9
1.25
*
1
0
-
06
1.35
*
1
0
-
03
1
/
128
4
7
1.84
*
1
0
-
08
3.23
*
1
0
-
04
1
/
256
3
5
2.82
*
1
0
-
10
7.99
*
1
0
-
05In this series of numerical experiments, both algorithms exhibit the same convergence properties as in the previous series of tests; that is, the number of iterations decreases as the number of grid points increases. From Table7, we can also see the sixth order convergence of the approximate solution to the exact solution of the boundary value problem on a sequence of grids. These results confirm the sixth order approximation of the compact explicit scheme proposed for the numerical implementation of the Sommerfeld-like boundary conditions.In Table8, we present the number of iterations required for the convergence of the GMRES and SKS methods for different coefficients k in the Helmholtz equation with the same boundary conditions used in the previous series of numerical tests. Data presented in Table 8 are similar to the data in Table 6.Table 8
3D Dirichlet-Sommerfeld problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
k
=
35.7
+
0.43
i
h
=
1
/
64
5
∣
7
6
∣
9
11
∣
43
16
∣
div
119
∣
div
14
∣
54
h
=
1
/
128
4
∣
5
4
∣
7
7
∣
8
7
∣
div
12
∣
11
7
∣
9
h
=
1
/
256
3
∣
4
3
∣
5
5
∣
6
5
∣
6
6
∣
7
5
∣
6
h
=
1
/
512
2
∣
4
3
∣
5
4
∣
5
3
∣
5
4
∣
6
4
∣
5
ψ
15
13
13
15
11
13In the last column of Table8, we present the results of calculations when k is a complex valued coefficient. The value of this coefficient approximately corresponds to the case of propagation of electromagnetic waves with the 1 GHz frequency in dry soil. These results confirm the effective implementation of the compact sixth order approximation scheme by using the proposed Krylov subspace type algorithms in the framework with the FFT-based low order preconditioners. The last row of the table suggests that the acceleration of the convergence with the increase of the number of grid points is much stronger than the acceleration in the previous series of experiments. It seems strongly dependent on the boundary conditions used in the numerical tests. In all our experiments, there is a clear connection between the parameter ψ and the improvement of convergence with the decrease of the step side. But the usefulness of this parameter in the analysis of the quality of a preconditioner requires further analysis since in majority of our numerical experiments the preconditioner does not satisfy the condition in Definition 1.The series of test problems considered suggests that in the majority of situations the preconditioned GMRES method is the most efficient choice for an effective implementation of the compact sixth order approximation scheme on the coarse grids but in the case of finer grids the SKS method in combination with lower-order approximation preconditioner presents an efficient alternative to the GMRES method. This alternative could become even more valuable when the GMRES method experiences stagnation.


## 4.1. 1D Test Problem

In the first series of numerical experiments, we consider the convergence of the compact sixth order approximation algorithm on a sequence of grids. In these tests we focus on the numerical solution of the 1D Dirichlet problem for the Helmholtz equation with zero boundary conditions on the interval0≤x≤1 and k=20. The computational grid in our experiments is defined by xi=ih,i=0,…,N+1, where h=1/(N+1) and N is the number of grid points. We consider the function u(xi)=xi(1-xi)cos⁡(kπxi),i=1,…,N as the exact solution of the original boundary value problem with the right-hand side f(xi)=-(2+k2(π2-xi(1-xi))cos⁡(kπxi) + 2kπ(2xi-1)sin(kπxi),i=1,…,N. The first column of Table 2 displays the step size h of the grids in our experiments. Columns 2–4 of the table present the number of iterations until the convergence required by the preconditioned GMRES method, the SKS method (Algorithm 1), and the Chebyshev acceleration (CA) algorithm [24]. In all our experiments we use restarted GMRES method with the restart value 20. We use the stoppage criterion ∥F-AU(n)∥2/∥F∥2≤10-10, where U(n) is the iterative solution on nth iteration. Column 5 reports the approximation error of the numerical solution Err6=max⁡1≤j≤N|Uj(n)-u(xj)|. In Column 6, the second order approximation error Err2 of the numerical solution obtained by using the preconditioner as a solver is presented for comparison. The number of grid points is chosen to satisfy the requirement h<2π/10k. The values of the minimum and maximum eigenvalues of the preconditioned matrix used by the Chebyshev acceleration algorithm are calculated exactly by (22).Table 2
1D Dirichlet problem.
Grid size (h)
GMRES
SKS
CA
Err
6
Err
2
1
/
32
6
8
6
3.88
*
1
0
-
03
1.16
*
1
0
-
01
1
/
64
4
4
5
5.08
*
1
0
-
05
2.4
*
1
0
-
02
1
/
128
3
3
4
7.61
*
1
0
-
7
5.7
*
1
0
-
03
1
/
258
3
3
3
1.14
*
1
0
-
8
1.4
*
1
0
-
03These experiments demonstrate the sixth order convergence of the numerical solution to the exact solution of the original boundary value problem. Theorem5 suggests that the residual on nth iteration in the presented iterative algorithms decreases with the increase in the number of grid points. The results of the numerical experiments confirm this conclusion. We have already shown that the order of the preconditioned matrix in this case is two. The numerical value of the constant ψ calculated by using (25) is 1.99. The 1D numerical experiments confirm the effectiveness of the iterative approaches under consideration and allow us to expect that the same properties hold for 3D problems.


## 4.2. 3D Test Problems

In the next series of experiments we consider different boundary value problems for the 3D Helmholtz equations (1) and (2). The computational domain is given by Ω={(x,y,z)∣0≤x,y,z≤1} with a uniform grid in all three directions xi=ih,i=1,…,Nx,yj=jh,j=1,…,Ny and zk=kh,k=1,…,Nz. The numbers of grid points Nx, Ny, and Nz in x,y, and z-directions are chosen to provide uniform grid size.Dirichlet Problem. First, we consider the Helmholtz equation with boundary conditions u=0 and k=20. This test is similar to the 1D case with Nx=Ny=Nz=N and h=1/(N+1). As the solution of the boundary value problem we choose the function u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy),ϕ(z)=sin(kπz), and 0≤x,y,z≤1.In the first series of tests we consider the convergence of the compact high order approximation finite-difference scheme (28) on a sequence of grids as well as convergence properties of the GMRES, SKS, and CA iterative algorithms. The iterative process is stopped when the initial residual is reduced by a factor of 10-10. The values of the minimum and maximum eigenvalues of the preconditioned matrix used by the Chebyshev acceleration algorithm are calculated exactly using (53). The description of the columns of Table 3 is similar to the description of Table 2.Table 3
3D Dirichlet problem.
Grid size (h)
GMRES
SKS
CA
Err
6
Err
2
1
/
64
5
10
23
3.03
*
1
0
-
04
4.47
*
1
0
-
06
1
/
128
3
6
15
8.13
*
1
0
-
05
6.35
*
1
0
-
08
1
/
256
3
5
14
2.06
*
1
0
-
05
9.68
*
1
0
-
10The test results confirm the sixth order approximation of the compact finite-difference scheme. Though, in our theoretical analysis (53) we could not prove the desired kth order preconditioning property in the case of the 3D Dirichlet problem, the SKS algorithm still exhibits acceleration of the convergence with the decrease of grid size. We also mention that the total computer time required for convergence for the SKS and GMRES methods is approximately the same and the Chebyshev acceleration algorithm requires about four times more computer time until convergence than the other two iterative algorithms do. The SKS algorithm has also greater potential for efficient implementation on parallel computers than the GMRES method does.One potential application of the proposed methods is the electromagnetic scattering problems. In these problems, the typical value of the coefficientk is in the range of 10 and 50. For example, the propagation of an electromagnetic wave with 1 GHz frequency in a vacuum is described by the Helmholtz equation with k≈20.9m-1. In Table 4 we present the convergence of the algorithms presented here on the sequence of grids with different values k. In every cell of the table we present the number of iteration for three methods under consideration: GMRES, SKS, and CA. The stoppage criterion is the same as in the previous series of experiments. If any of these three methods fails to converge in 100 iterations, we indicate this by using symbol “>100" and in the case of divergence we use “div” symbol.Table 4
3D Dirichlet problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
4
∣
7
∣
14
5
∣
10
∣
23
6
∣
13
∣
25
12
∣
div
∣ >100
19
∣
div
∣ >100
h
=
1
/
128
3
∣
5
∣
14
3
∣
6
∣
15
4
∣
8
∣
17
4
∣
10
∣
69
5
∣
11
∣ >100
h
=
1
/
256
3
∣
4
∣
12
3
∣
5
∣
14
3
∣
6
∣
13
3
∣
6
∣
22
4
∣
7
∣
63
ψ
2.02
2.00
1.98
2.00
2.02From these numerical experiments we can see that the preconditioned GMRES and SKS methods demonstrate excellent convergence properties. We also observe that the convergence of the proposed algorithms improves with the increase in the number of the grid points. In the analysis of the algorithm’s convergence (54), it was shown that the number of iterations does not increase with the increase of the grid points. But the last row in Table 4 indicates that the parameter ψ (25) for all values of k is approximately 2, which indicates that similar to the 1D case the convergence of the 3D algorithm not only does not depend on the grid size but accelerates with the increase of the number of grid points. This also indicates that even in the case when the SKS method does not converge, one can expect convergence on a finer grid. In all numerical tests, the parameter ψ is calculated by using the first two iterations of the SKS method on the grids with the grid sizes h=1/128 and h=1/256.From the results presented, it is clear that slow convergence and difficulties arising in the choice of the parametersm^ and M^ in the Chebyshev acceleration algorithm (18) make this method a poor choice for the implementation of the developed high-order approximation scheme. These numerical tests illustrate the fact that the Chebyshev polynomial is not optimal on the discrete spectrum and this has a dramatic effect on the iterative methods based on this polynomial. So, in the next numerical experiments we focus on the convergence properties of the GMRES and SKS methods in our numerical framework.Dirichlet-Neumann Boundary Conditions. In the next series of experiments we replace the Dirichlet boundary conditions u=0 in (2) with u=ub at the top of the computational domain (z=1) and with the Neumann boundary conditions ∂u/∂z=0 at z=0. On the other boundaries of the rectangular computational domain Ω={0≤x,y,z≤1} we still use the Dirichlet boundary conditions u=0. As a test function that satisfies the boundary conditions we consider u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy),ϕ(z)=cos⁡(kπz), and 0≤x,y,z≤1. Also, ub=ϕ1(x)ϕ2(y)ϕ3(1),0≤x,y≤1. The number of grid points is chosen to provide the same grid step in all three directions; that is Nx=Ny=Nz-1. The stoppage criterion is the same as in the previous numerical experiments. Table 5 provides the results of the numerical tests on a sequence of grids with the coefficient k=20. As already mentioned we focus only on the two best algorithms: GMRES and SKS. The rest of Table 5 is similar to Table 3.Table 5
3D Dirichlet-Neumann problem.
Grid size (h)
GMRES
SKS
Err
6
Err
2
1
/
64
8
11
4.65
*
1
0
-
06
3.67
*
1
0
-
04
1
/
128
5
8
6.61
*
1
0
-
08
8.74
*
1
0
-
05
1
/
256
4
7
1.01
*
1
0
-
09
2.17
*
1
0
-
05The results of this series of numerical tests indicate the sixth order convergence of the resulting compact finite-difference scheme including explicit approximation of the Neumann boundary conditions. As expected, the number of iterations decreases with the increase of the number of grid points in both the GMRES and SKS approaches. In Table6 we present the number of iterations required for convergence for both algorithms on a sequence of grids and with different k coefficients in the Helmholtz equation. The content of Table 6 is similar to the content of Table 4, except for the number of iterative approaches under consideration. We present only the number of iterations for the GMRES and SKS methods.Table 6
3D Dirichlet-Neumann problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
h
=
1
/
64
6
∣
10
8
∣
11
18
∣
div
36
∣
div
>600∣div
h
=
1
/
128
4
∣
9
5
∣
8
8
∣
9
14
∣
div
28
∣
div
h
=
1
/
256
3
∣
8
4
∣
7
5
∣
div
9
∣
div
7
∣
div
h
=
1
/
512
2
∣
6
3
∣
6
3
∣
6
5
∣
6
4
∣
6
ψ
2.57
2.18
2.39
1.65
1.94It follows from this table that the GMRES approach is much more efficient on coarse grids but on finer grids, the processor time for the SKS method is the same or even smaller than for the GMRES algorithm. For example, on the grid withh=1/512 and k=50 the corresponding processor times for the GMRES and SKS methods are 1467 sec and 1052 sec. Also, as we mentioned before, the SKS method has much greater potential than the GMRES algorithm for efficient implementation on parallel computers. We also expect that even if the SKS method does not converge on coarser grids, by reducing the grid size one can achieve convergence of this algorithm. We can see this in the example of the convergence of the methods under consideration in the case of k≥30. The values of the parameter ψ vary significantly for different values of k but ψ>0 still indicates that the convergence of the iterative algorithms improves with the increase in the number of the grid points.Dirichlet-Sommerfeld-Like Boundary Conditions. In the last series of numerical experiments, we consider the boundary value problem for the Helmholtz equation (1) with a combination of Sommerfeld-like (radiation) boundary conditions (2) at z=0 and z=1, and Dirichlet boundary conditions at all other boundaries of the rectangular computational domain Ω={0≤x,y,z≤1}. The source function f in (1) selected that the true solution is u(x,y,z)=ϕ1(x)ϕ2(y)ϕ3(z), where ϕ1(x)=x3(1-x)3,ϕ2(y)=y(1-y)cos⁡(kπy), ϕ3(z)=eikz+e-ik(z-1)-2, and 0≤x,y,z≤1. Note that the analytic solution satisfies the radiation boundary condition (2). The numbers of grid points in x,y,z-directions are chosen such that Nx=Ny=Nz-2. First we consider the convergence of our sixth order approximation scheme on the sequence of grids with k=20. The main goal of this series of numerical experiments is to investigate the convergence properties of the new explicit compact sixth order approximation of the Sommerfeld-like boundary conditions (47) proposed in this paper. As in previous test runs, the iterative process is stopped when the l2-norm of the initial residual is reduced by a factor of 10-10. The results are presented in Table 7. The description of data presented in this table is the same as in the case of Table 5.Table 7
3D Dirichlet-Sommerfeld problem.
Grid size (h)
GMRES
SKS
Err
6
Err
2
1
/
64
6
9
1.25
*
1
0
-
06
1.35
*
1
0
-
03
1
/
128
4
7
1.84
*
1
0
-
08
3.23
*
1
0
-
04
1
/
256
3
5
2.82
*
1
0
-
10
7.99
*
1
0
-
05In this series of numerical experiments, both algorithms exhibit the same convergence properties as in the previous series of tests; that is, the number of iterations decreases as the number of grid points increases. From Table7, we can also see the sixth order convergence of the approximate solution to the exact solution of the boundary value problem on a sequence of grids. These results confirm the sixth order approximation of the compact explicit scheme proposed for the numerical implementation of the Sommerfeld-like boundary conditions.In Table8, we present the number of iterations required for the convergence of the GMRES and SKS methods for different coefficients k in the Helmholtz equation with the same boundary conditions used in the previous series of numerical tests. Data presented in Table 8 are similar to the data in Table 6.Table 8
3D Dirichlet-Sommerfeld problem (differentk2).
k
=
10
k
=
20
k
=
30
k
=
40
k
=
50
k
=
35.7
+
0.43
i
h
=
1
/
64
5
∣
7
6
∣
9
11
∣
43
16
∣
div
119
∣
div
14
∣
54
h
=
1
/
128
4
∣
5
4
∣
7
7
∣
8
7
∣
div
12
∣
11
7
∣
9
h
=
1
/
256
3
∣
4
3
∣
5
5
∣
6
5
∣
6
6
∣
7
5
∣
6
h
=
1
/
512
2
∣
4
3
∣
5
4
∣
5
3
∣
5
4
∣
6
4
∣
5
ψ
15
13
13
15
11
13In the last column of Table8, we present the results of calculations when k is a complex valued coefficient. The value of this coefficient approximately corresponds to the case of propagation of electromagnetic waves with the 1 GHz frequency in dry soil. These results confirm the effective implementation of the compact sixth order approximation scheme by using the proposed Krylov subspace type algorithms in the framework with the FFT-based low order preconditioners. The last row of the table suggests that the acceleration of the convergence with the increase of the number of grid points is much stronger than the acceleration in the previous series of experiments. It seems strongly dependent on the boundary conditions used in the numerical tests. In all our experiments, there is a clear connection between the parameter ψ and the improvement of convergence with the decrease of the step side. But the usefulness of this parameter in the analysis of the quality of a preconditioner requires further analysis since in majority of our numerical experiments the preconditioner does not satisfy the condition in Definition 1.The series of test problems considered suggests that in the majority of situations the preconditioned GMRES method is the most efficient choice for an effective implementation of the compact sixth order approximation scheme on the coarse grids but in the case of finer grids the SKS method in combination with lower-order approximation preconditioner presents an efficient alternative to the GMRES method. This alternative could become even more valuable when the GMRES method experiences stagnation.


## 5.  Conclusions

New 3D compact sixth order explicit finite-difference schemes for the approximation of Neumann and Sommerfeld-like boundary conditions on rectangular computational domains with uniform grid size were developed and implemented. Together with the compact sixth order approximation scheme for the Helmholtz equation proposed in [6], these algorithms represent highly accurate methods for the solution of boundary value problems for Helmholtz equations.A new rapid iterative method based on preconditioned Krylov subspace methodology was developed for the implementation of the proposed compact finite-difference schemes. The strategy is based on a combination of higher order approximation schemes and a lower-order approximation preconditioner. The analysis of some typical test problems reveals the attractive properties of the developed methods such as the decrease of the number of iterations until convergence with the increase of the number of the grid points or the size of the resulting matrix. This approach is especially attractive in situations in which the lower approximation solver already exists and the original boundary value problem calls for more accurate approximation.The typical time to produce the sixth order accuracy solution of the 3D Helmholtz equation with a combination of the Dirichlet and Sommerfeld boundary conditions on a5123 grid was just 30 minutes on an iMac using only one 2.93 GHz Intel Core i7 processor. The method was tested for realistic parameter ranges typical for electromagnetic scattering problems. We must notice that the Sommerfeld-like boundary is just a first order approximation of the Sommerfeld conditions on the unbounded domain. So, direct application of the higher order approximation for the Sommerfeld-like boundary condition to the solution of a scattering problem is not always justified. However, a straightforward extension of the approximation approach presented in this paper could be applied for approximation of the absorbing boundary conditions (see, e.g., [13]) or can be used in the implementation of the perfectly matched layer (PML) boundary conditions (see, e.g., [25]).



---

*Source: 745849-2014-01-21.xml*