# A Robust Optimal Finite Difference Scheme for the Three-Dimensional Helmholtz Equation

Authors: Dongsheng Cheng; Baowen Chen; Xiangling Chen
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

We propose a robust optimal 27-point finite difference scheme for the Helmholtz equation in three-dimensional domain. In each direction, a special central difference scheme with 27 grid points is developed to approximate the second derivative operator. The 27 grid points are divided into four groups, and each group is involved in the difference scheme by the manner of weighted combination. As for the approximation of the zeroth-order term, we use the weighted average of all the 27 points, which are also divided into four groups. Finally, we obtain the optimal weights by minimizing the numerical dispersion with the least-square method. In comparison with the rotated difference scheme based on a staggered-grid method, the new scheme is simpler, more practical, and much more robust. It works efficiently even if the step sizes along different directions are not equal. However, rotated scheme fails in this situation. We also present the convergence analysis and dispersion analysis. Numerical examples demonstrate the effectiveness of the proposed scheme.


---

## Body

## 1. Introduction

Many of the physical problems are governed by the wave equation, which is also well known as the Helmholtz equation in frequency domain. The Helmholtz equation finds a wide application in many fields of science, engineering, and industry. To solve the Helmholtz equation numerically, finite difference methods [1–3]) and finite element methods [4–7] are frequently employed. For its easy implementation and less computational complexity, the finite difference method is usually preferred, especially in the engineering field such as oil-gas exploration.All of the numerical methods suffer from the so-called “pollution effect” [8], which leads to the dispersion (phase) errors. In the cases of high-dimensional problems, the “pollution effect” cannot be eliminated; then, the dispersion errors always exist. Consequently, a good numerical solver for the Helmholtz equation should have a good performance in reducing the dispersion errors. For finite difference methods, the classical central difference scheme gives a poor performance with producing large errors. To improve the performance, Jo et al. [9] proposed the rotated 9-point difference scheme, based on the staggered-grid method [10]. The rotated 9-point scheme combined the classical Cartesian coordinate system and its rotated systems to discretize the two-dimensional (2D) Helmholtz equation; that is, the discretization of the equation is a weighted combination of discretization on the Cartesian system and that on the rotated system. The optimal weights were determined by minimizing the dispersion errors. The rotated 9-point scheme obtained a good effect and was extended to other situations (cf. [11–13]). In [11], Shin et al. constructed a 2-D 25-point rotated difference scheme with determining the optimal weights by the singular-value decomposition method. In [12], Hustedt et al. developed a group of new finite difference methods based on staggered-grid stencils of various accuracy and grid rotation strategies. For the three-dimensional (3-D) Helmholtz Equation, [13] employed staggered-grid stencils on 7 rotated coordinate systems to approximate the differential operators, which is rather complicated. Motivated by the rotated schemes, many weighted-averaging difference schemes were designed by combining the central second-order difference scheme at the central node and its neighboring nodes (cf. [14–16]). High order schemes [1, 17, 18] are also established to improve the accuracy; however, they impose strict requirements on the smoothness of the right-hand side.The rotated difference schemes are not robust in the sense that they fail when the step sizes are not equal in different directions. Hence, the rotated schemes cannot handle the problems with nonequidistant sampling in different directions. Furthermore, the rotated schemes are complex, especially for the 3D Helmholtz equation, which is frequently used in practical applications. To obtain a robust difference scheme for the 3D Helmholtz equation, in this paper, we propose a robust optimal 27-point difference scheme, which remains weighted and second order in accuracy but is rotation-free. To discretize the second derivative operator of the 3D Helmholtz equation, the new scheme utilizes a generalized central difference formula, which is obtained by replacing each of the 7 grid points with a weighted combination of it and its 11 neighboring grid points. To discretize zeroth-order term, we employ the weighted average of the 27 points. Both the discretizations of the second derivative operator and the zeroth-order term use all of the 27 grid points, which are split into four groups. To determine the weights, the least-square method is used to solve the optimization problem of minimizing the dispersion error. In comparison with the rotated scheme, the new scheme is simple, easy to use. More importantly, it is robust enough to handle the problems with nonequidistant sampling, while the rotated scheme fails in this situation. We prove this by theoretical analysis as well as numerical experiments.We point out that both the new scheme and rotated scheme are similar to the “dispersion relation preserving” (DRP) schemes that were originally introduced in [19] and continued with [20], since all of these methods are based on reducing the numerical dispersion error related to wave problems. However, the DRP schemes are usually suitable for the finite-difference time-domain (FDTD) simulations. For instance, [21] developed the DRP algorithms to minimize the numerical dispersion error in large-scale three-dimensional FDTD simulations. In [22], the DRP method was further generalized to include specification of both dispersion behavior and the orders of accuracy (OoAs) for both the wave equation (WE) and Maxwell’s equations (MEs). Based on DRP techniques, [23] presented higher order FDTD simulations with controllable dispersion error. In [24], a general procedure was introduced to reduce grid dispersion error in broadband FDTD simulations. In [25], a unified methodology for FDTD simulation was presented based on certain modifications of the characteristic equation that accompanies any given discretized version of the WE. In [26], this development was expanded into a comprehensive new methodology for the systematic generation of WE-FDTD schemes tailored to the spectrum of the excitation. In this paper, the proposed method differs from the DRP schemes since it is formulated to solve the Helmholtz equation only in spatial domain.The remainder of this paper is organized as follows. In Section2, we present a robust optimal 27-point finite difference scheme for the 3D Helmholtz Equation and prove it to be second order in accuracy. In Section 3, we perform dispersion error analysis and obtain a group of optimal weights for the new scheme. In Section 4, numerical experiments are present to validate the efficiency and effectiveness of the new scheme. Finally, Section 5 concludes this paper.


## 2. A Robust Optimal 27-Point Finite Difference Scheme

In this section, we propose a new finite difference scheme for the 3D Helmholtz equation. We also present the convergence analysis, which shows that new scheme is second order in accuracy.The 3D Helmholtz equation reads(1)-∂2u∂x2-∂2u∂y2-∂2u∂z2-k2u=ginR3,where u denotes the pressure wavefield, k≔2πf/v is the wavenumber with f, v are the frequency and velocity, respectively, and g denotes a source function.To describe the new finite difference scheme, we use the network of grid points(xm,yn,zl). Here, xm≔x0+(m-1)hx, yn≔y0+(n-1)hy, and zl≔z0+(l-1)hz, with hx,hy,hz being the step sizes in the x-, y-, and z-directions, respectively. Note that, for the situation of equidistant sampling in all three directions, we have hx=hy=hz, while for nonequidistant sampling, not all of them are equal. We next present the 27-point finite difference stencil with numbering in Figure 1, where the grid points are numbered as (m0,n0,l0) with m0,n0,l0∈{-1,0,1}. Moreover, (0,0,0) represents the center point of the stencil, and the others denote the points surrounding it. For simplicity, the discretization of a function φ at point (m0,n0,l0) is denoted by φm0n0l0 or φm0n0l0.Figure 1
The 27-point finite difference stencil with numbering.The discretization of the 3D Helmholtz equation comprises two segments, one of which is the discretization of the second derivative operator∂2u/∂x2+∂2u/∂y2+∂2u/∂z2 and the other is the discretization of the zeroth-order term k2u. To approximate ∂2u/∂x2,∂2u/∂y2 and ∂2u/∂z2, we establish a generalized central difference scheme as follows:(2)Lxu000≔1hx2u100R-u0,0,0R-u000L-u-100L,(3)Lyu000≔1hy2u010F-u0,0,0F-u000B-u0-10B,(4)Lzu000≔1hz2u001D-u000D-u000U-u00-1U,where (5)u100R≔a1∑l∈Ωu10l+∑n∈Ωu1n0+u000+a2∑n,l∈Ωu1nl+∑l∈Ωu00l+a3u100,u000R≔a1∑l∈Ωu00l+∑n∈Ωu0n0+u100+a2∑n,l∈Ωu0nl+∑l∈Ωu10l+a3u000,u-100L≔a1∑l∈Ωu-10l+∑n∈Ωu-1n0+u000+a2∑n,l∈Ωu-1nl+∑l∈Ωu-10l+a3u-100,u000L≔a1∑l∈Ωu00l+∑n∈Ωu0n0+u-100+a2∑n,l∈Ωu0nl+∑l∈Ωu-10l+a3u000,u010F≔a1∑l∈Ωu01l+∑m∈Ωum10+u000+a2∑m,l∈Ωum1l+∑l∈Ωu00l+a3u010,u000F≔a1∑l∈Ωu00l+∑m∈Ωum00+u010+a2∑m,l∈Ωum0l+∑l∈Ωu01l+a3u000,u0-10B≔a1∑l∈Ωu0-1l+∑m∈Ωum-10+u000+a2∑m,l∈Ωum-1l+∑l∈Ωu00l+a3u0-10,u000B≔a1∑l∈Ωu00l+∑m∈Ωum00+u0-10+a2∑m,l∈Ωum0l+∑l∈Ωu0-1l+a3u000,u001D≔a1∑n∈Ωu0n1+∑m∈Ωum01+u000+a2∑m,n∈Ωumn1+∑n∈Ωu0n0+a3u001,u000D≔a1∑n∈Ωu0n0+∑m∈Ωum00+u001+a2∑m,n∈Ωumn0+∑n∈Ωu0n1+a3u000,u00-1U≔a1∑n∈Ωu0n-1+∑m∈Ωum0-1+u000+a2∑m,n∈Ωumn-1+∑n∈Ωu0n0+a3u00-1,u000U≔a1∑n∈Ωu0n0+∑m∈Ωum00+u00-1+a2∑m,n∈Ωumn0+∑n∈Ωu0n-1+a3u001,with Ω≔{1,-1} and parameters a1,a2,a3 being the weights. As is observed, (2)-(4) reserve the form of classical central difference scheme; however, instead of using a single grid-point, a set of grid points are employed. The grid points consist of a point and its 11 neighbors in the manner of weighted combination; for instance, u1,0,0R is a combination of u1,0,0 and u101, u10-1, u110, u1-10, u000, u111, u11-1, u1-11, u1-1-1, u001, u00-1. In the x-direction, all the 27 points are divided into four groups which correspond to u100R, u000R, u-100L, and u000L, respectively. For the y- and z-directions, the situations are similar.To approximate the term of zero orderk2u, we utilize the weighted average of all the 27 points as follows:(6)Ik2u000≔b1k2u000+b26∑m∈Ωk2um00+∑n∈Ωk2u0n0+∑l∈Ωk2u00l+b312∑m,n∈Ωk2umn0+∑n,l∈Ωk2u0nl+∑m,l∈Ωk2um0l+b48∑m,n,l∈Ωk2umnl,where parameters b1,b2,b3,b4 are another group of weights. As is seen, all the 27 points are involved in the discretization of k2u, and they are also divided into four groups, which correspond to b1,b2,b3,b4 respectively.Finally, the discretization of equation (1) at point (0,0,0) is expressed as(7)-Lxu000-Lyu000-Lyu000-Ik2u000=g000,which gives a new 27-point finite difference scheme for the 3D Helmholtz equation. Substituting (2), (3), (4), and (6) into (7), the new scheme is rewritten as(8)∑m,n,l∈ΦCmnlumnl=g000,where Φ≔{-1,0,1} and Cmnl are the coefficients with(9)Cmnl=T1≔-a2hx+a2hy+a2hz-b48k2,form,n,l∈ΩT2≔2a2hx+a2-a3hy+a2-a1hz-b312k2,form=0,n∈Ω,l∈ΩT3≔a2-a1hx+2a2hy-a1hz-b312k2,form∈Ω,n=0,l∈ΩT4≔2a2-a1hx+2a1-a2hy+a1-a3hz-b26k2,form=0,n=0,l∈ΩT5≔-a1hx+a1hy-2a2hz-b312k2,form∈Ω,n∈Ω,l=0T6≔2a1hx+a1-a3hy+2a1-a2hz-b26k2,form=0,n∈Ω,l=0∈ΩT7≔a1-a3hx+2a1hy+2a1hz-b26k2,form∈Ω,n=0∈Ω,l=0T8≔2a3-a1hx+2a3-a1hy+2a3-a1hz-b1k2,form=0,n=0,l=0.We next present the convergence analysis of the new difference scheme. For the new scheme (7), we have the following proposition.Proposition 1.
The new difference scheme (7) is second order in accuracy if the weights a1,a2,a3,b1,b2,b3,b4 satisfy(10)3a1+2a2+a3=1,b1+b2+b3+b4=1.Proof.
We start with handling the formula (2). Let P1≔u100R-u000R, P2=u-100L-u000L, r1≔a1-a2, r2≔a1-a3; then, we have(11)P1=r1∑l∈Ωu10l-u00l+a1∑n∈Ωu1n0-u0n0+a2∑n,l∈Ωu1nl-u0nl-r2u100-u000,(12)P2=r1∑l∈Ωu-10l-u00l+a1∑n∈Ωu-1n0-u0n0+a2∑n,l∈Ωu-1nl-u0nl)-r2u-100-u000.Let hx=h,hy=αh,hz=βh. Applying the Taylor expansion to each term of (11) and (12) yields(13)P1=γu100-u000+h3a1∂B∂x1/200+2αβa2∂A∂x1/200-a2∂∂x∂2u∂z21/200+oh4,(14)P1=γu100-u000+h3a1∂B∂x-1/200+2αβa2∂A∂x-1/200-a2∂∂x∂2u∂z2-1/200+oh4,where γ≔3a1+2a2+a3, A≔∂2u/∂y2+∂2u/∂z2, and B≔α2∂2u/∂y2+β2∂2u/∂z2. Substituting (14) and (13) into (2) and applying the Taylor expansion again, we obtain(15)Lxu000=γ∂2u∂x2000+h2D+oh2,where(16)D≔124∂4u∂x4000+2αβa2+α2a1∂4u∂x2∂y2000+2αβa2+β2a1-β2a2∂4u∂x2∂z2000.
Similarly, for formulas (3) and (4), we have(17)Lyu000=c∂2u∂y2000+h2E+oh2.(18)Lzu000=c∂2u∂z2000+h2F+oh2.
where(19)E≔α424∂4u∂y4000+2βa2+αa1∂4u∂y2∂z20002βa2+αβ2a1-αβ2a2∂4u∂y2∂x2000,F≔β224∂4u∂z4000+2αa2+βa1∂4u∂z2∂x2000+2αa2+βα2a1-α2a2∂4u∂z2∂y2000.
For the formula (6), denote its second, third, and fourth parts by Q1,Q2, and Q3 respectively. According to the Taylor expansion, we obtain(20)Q1=b2k2u000+b26h2∂2∂x2k2u000+α2∂2∂y2k2u000+β2∂2∂z2k2u000+oh3,(21)Q2=b3u000+b33h2∂2∂x2k2u000+α2∂2∂y2k2u000+β2∂2∂z2k2u000+oh3,(22)Q3=b4u000+b48h2∂2∂x2k2u000+2α2∂2∂y2k2u000+4β2∂2∂z2k2u000+oh3,Finally, substituting (15), (17), (18), (20), (21), and (22) into the left-hand side of the difference scheme (7) leads to(23)-Lxu000-Lyu000-Lyu000-Ik2u000=3a1+2a2+a3∂2u∂x2+∂2u∂y2+∂2u∂z2000+b1+b2+b3+b4k2u000+h2D+E+F+H+oh2,with(24)H≔b26+b33+b48∂2∂x2k2u000+α2b26+b33+b44∂2∂y2k2u000+β2b26+b33+b42∂2∂z2k2u000.
As is observed from (23), the new difference scheme has a second-order convergence rate when 3a1+2a2+a3=1 and b1+b2+b3+b4=1.Assume3a1+2a2+a3=1, b1+b2+b3+b4=1. Then, according to Proposition 1, the new difference scheme is second order in accuracy. In the next section, we present the determination of the weights a1,a2,a3,b1,b2,b3,b4 based on minimizing the dispersion error.


## 3. Dispersion Error Analysis and Determination of the Weights for the New Difference Scheme

### 3.1. Dispersion Error Analysis

In this subsection, we first perform dispersion analysis for the new difference scheme, which reveals the error between the exact wavenumberk and the numerical wavenumber obtained from the new scheme. The dispersion error indicates that the new scheme is a robust one, which can handle the problems of nonequidistant sampling in different dimensions.Assumehx=h, hy=αh, hz=βh, which gives flexible step sizes in different dimension. When α=β=1, it corresponds to equidistant sampling; otherwise, it corresponds to nonequidistant sampling. Moreover, let 3a1+2a2+a3=1, b1+b2+b3+b4=1, and k be a constant wavenumber.For the 3D Helmholtz equation (1) with g=0, its exact solution is(25)ux,y,z≔eikxcos⁡ϕcos⁡θ+ycos⁡ϕsin⁡θ+zsin⁡ϕ,which is also known as the classical plane-wave solution. Here, i is the imaginary unit and θ and ϕ are the propagation angles from the x-axis and z-axis, respectively. On the difference stencil, the discrete form of (25) is given by(26)umnl≔eikhcos⁡ϕcos⁡θ+αcos⁡ϕsin⁡θ+βsin⁡ϕ.Let λ denote the wavelength; then, λ/hx, λ/hy and λ/hz are the numbers of grid points per wavelength in the directions of x-, y-, and z-axis, respectively. For convenience, choose G≔λ/hx=λ/h as a benchmark. It is well known that λ=v/f. Together with k=2πf/v, we have kh=2π/G. Then, substituting the discrete plane-wave solution umnl into difference scheme (8) with g000=0, we obtain(27)8abcT1+4bcT2+4acT3+2cT4+4abT5+2bT6+2aT7+T8=0,where(28)a≔cos2πGcos⁡ϕcos⁡θ,b≔cosα2πGcos⁡ϕsin⁡θ,c≔cosβ2πGsin⁡ϕ.Let kN denote the numerical wavenumber, and replace k with kN in Tj (j=1,2,…,8). Then, some straightforward manipulations yield(29)kNk=G2πMN,where(30)N≔b1+b23a+b+c+b33ab+bc+ac+abcb4,M≔-8abca21+s+t+4bc[2a2-r1s+t+4ac2sa2-a1t-r1-4aba1+a1s-2a2t+2c2r11+s+r2t+2b2a1+r2s+2r1t+2ar2+2a1s+2a1t-2r21+s+t, with s=1/α2,t=1/β2, r1≔a1-a2, r2≔a1-a3.Definep1≔cos⁡ϕ,p2≔sin⁡ϕ,q1≔cos⁡θ,q2≔sin⁡θ, and (31)L1≔-6α2p12q22-6β2p22+2p12q12-8α2+β2p1p2-8p1q1,L2≔-40β2p1p2-32(1+sβ2p1p2q1-481+α2p12q1q2-8α24β2+5α2p1p2q2-4β2p22-4p12q12-4α2p12q22-48tα2p12q1q2,L3≔-2β2p22-2α2p12q22-2p12q12,L4≔b2+2b3+3b43p12q12+α2p12q22+β2p22.Then, we have the following proposition which reveals the error between the numerical wavenumber kN and the exact wavenumber k.Proposition 2.
For the 27-point finite difference scheme (7), the numerical wavenumber kN and the exact wavenumber k satisfy(32)kN=k+12a1L1+a2L2+a3L3+L4k3h2+Ok4h3,kh→0.Proof.
Letτ≔kh, and rewrite a,b,c as(33)aτ=cosτcos⁡ϕcos⁡θ,b=cosατcos⁡ϕsin⁡θ,c=cosβτsin⁡ϕ.Hence, M,N can also be considered as two functions with respect to τ. Then, we compute the Taylor expansions of M(τ),1/N(τ) around 0, and after some straightforward manipulations, we obtain(34)Mτ=τ2+τ4a1L1+a2L2+a3L3+Oτ5,τ→0(35)1Nτ=1+τ2L4+Oτ3.τ→0Inserting (34), (35) into (29) yields(36)kNk=1+τ2a1L1+a2L2+a3L3+L4+Oτ3,τ→0Finally, computing again the Taylor expansion of the right-hand side of (36) and replacing τ with kh, we obtain formula (32), which completes the proof.As can be seen from Proposition2, kN converges to k in a second order, which is independent of parameters α and β. This is important since it means that the convergence is robust even if the step sizes in three dimensions are completely different. However, for the 3D rotated difference scheme, it holds(37)kN=μα,βk+νk3h2+Ok4h3,kh→0.where μ(α,β)=1 only for α=β=1 (equidistant sampling). When α≠1 or β≠1 (nonequidistant sampling), μ(α,β)≠1, which means kN does not converge to k, and the rotated scheme fails. Hence, Proposition 2 shows that the new scheme is robust, namely, work effectively for both the situations of equidistant and nonequidistant sampling in different dimensions.


### 3.2. Determination of the Weights

In this subsection, we determine the weights for the new scheme based on minimizing the dispersion error, which is transformed to an optimization problem. We use the least-square method to solve the optimization problem, obtaining a group of optimal weights.For formula (29) in Section 3.1, it is also known as the dispersion relation formula, which is equivalent to the normalized phase velocity [27, 28]. Ideally, kN/k equals 1, which means that the difference scheme is exact. However, for the 2D and 3D situations, it is impossible to achieve this due to the “pollution effect”. Hence, we expect kN/k could be close to 1. For this end, set(38)Ja1,a2,a3,b1,b2,b3,b4;G,ϕ,θ≔kNk-1.When determining the weight parameters a1,a2,a3,b1,b2,b3,b4, (38) is expected to be close to 0. This can be transformed to the following optimization problem:(39)minJa1,a2,a3,b1,b2,b3,b4;G,ϕ,θ∞,s.t.3a1+2a2+a3=1,∑j=14bj=1,Gmin≤G≤Gmax,0≤ϕ≤π2,0≤θ≤π2,where Gmax is usually set to be 400 and Gmin is not smaller than 2 according to the Nyquist sampling limit. To handle (39), let J(a1,a2,a3,b1,b2,b3,b4;G,ϕ,θ)=0, and then obtain(40)A1a1+A2a2+B1b1+B2b2+B3b3=S,where (41)A1≔4bcs+t+4ac1+t+4ab1+s-4a1+r-4bs+r-4ct+r+8r,A2≔8abcr-8abt-4ac1+2s-4bc1+r-4a-4bcs-t-4ct-s-1+4r,B1≔4π2G21-abc,B2≔4π2G213a+b+c-abc,B3≔4π2G213ab+bc+ac-abc,S≔2r-2a+bs+ct-abc4π2G2,with s=1/α2,t=1/β2, r≔1+s+t. Next, we take equidistant samplings for variables G,ϕ and θ in the corresponding intervals Gmin,Gmax,[0,π/2] and [0,π/2] respectively. Assume the numbers of samples for G,ϕ and θ are m,n, and l; then, inserting the samples into (40) yields an overdetermined linear system with unknowns being a1,a2,b1,b2,b3 and size being m×n×l. We use the least-square method to solve the overdetermined linear system and finally obtain a group of optimal weights. For different α and β, the weights vary accordingly. It means that the new difference scheme has a unified form, but equipped with flexible weights, which adaptively vary with the step sizes in three dimensions. In Table 1, we present several groups of weights corresponding to different α and β.Table 1
Weights corresponding to differentα and β.
( α ,  β ) a 1 a 2 a 3 b 1 b 2 b 3 b 4 (1.0, 1.0) 0.1009 0.0011 0.6952 0.3714 0.8111 -0.2928 0.1102 (1.0, 1.5) 0.1154 0.0112 0.6313 0.2614 0.9925 -0.4020 0.1481 (1.5, 1.0) 0.1001 0.0007 0.6984 0.4242 0.6618 -0.1515 0.0656 (1.5, 1.5) 0.1194 0.0061 0.6297 0.3340 0.7907 -0.2159 0.0912 (1.0, 2.0) 0.1522 0.0247 0.4939 0.0274 1.3965 -0.6477 0.2238 (2.0, 1.0) 0.1029 -0.0016 0.6945 0.4275 0.6420 -0.1202 0.0508 (2.0, 2.0) 0.1619 0.0051 0.5040 0.2015 0.9602 -0.2963 0.1346The penalty function in optimization problem (39) is related to the difference between the exact and numerical wavenumbers. Since penalizing all wavenumbers equally, this may not be the most appropriate one. In practice, the more appropriate would be one that penalizes more strongly the error at large wavenumbers, since the cumulative dispersion error is more critical, and the spectral pollution problem is more acute for higher frequencies/wavenumbers. To remedy this, we adopt a flexible strategy for determination of the weights; that is, we estimate the interval IG≔[Gmin,Gmax] by using a priori information. In practice, for example, if the frequency f∈[fmin,fmax] and the velocity v∈[vmin,vmax], then, for a given step size h, we have Gmin≔vmin/hfmax and Gmax≔vmax/hfmin. Hence, the interval IG is flexible according to different cases, and a group of appropriate weights can be obtained for the difference scheme.


## 3.1. Dispersion Error Analysis

In this subsection, we first perform dispersion analysis for the new difference scheme, which reveals the error between the exact wavenumberk and the numerical wavenumber obtained from the new scheme. The dispersion error indicates that the new scheme is a robust one, which can handle the problems of nonequidistant sampling in different dimensions.Assumehx=h, hy=αh, hz=βh, which gives flexible step sizes in different dimension. When α=β=1, it corresponds to equidistant sampling; otherwise, it corresponds to nonequidistant sampling. Moreover, let 3a1+2a2+a3=1, b1+b2+b3+b4=1, and k be a constant wavenumber.For the 3D Helmholtz equation (1) with g=0, its exact solution is(25)ux,y,z≔eikxcos⁡ϕcos⁡θ+ycos⁡ϕsin⁡θ+zsin⁡ϕ,which is also known as the classical plane-wave solution. Here, i is the imaginary unit and θ and ϕ are the propagation angles from the x-axis and z-axis, respectively. On the difference stencil, the discrete form of (25) is given by(26)umnl≔eikhcos⁡ϕcos⁡θ+αcos⁡ϕsin⁡θ+βsin⁡ϕ.Let λ denote the wavelength; then, λ/hx, λ/hy and λ/hz are the numbers of grid points per wavelength in the directions of x-, y-, and z-axis, respectively. For convenience, choose G≔λ/hx=λ/h as a benchmark. It is well known that λ=v/f. Together with k=2πf/v, we have kh=2π/G. Then, substituting the discrete plane-wave solution umnl into difference scheme (8) with g000=0, we obtain(27)8abcT1+4bcT2+4acT3+2cT4+4abT5+2bT6+2aT7+T8=0,where(28)a≔cos2πGcos⁡ϕcos⁡θ,b≔cosα2πGcos⁡ϕsin⁡θ,c≔cosβ2πGsin⁡ϕ.Let kN denote the numerical wavenumber, and replace k with kN in Tj (j=1,2,…,8). Then, some straightforward manipulations yield(29)kNk=G2πMN,where(30)N≔b1+b23a+b+c+b33ab+bc+ac+abcb4,M≔-8abca21+s+t+4bc[2a2-r1s+t+4ac2sa2-a1t-r1-4aba1+a1s-2a2t+2c2r11+s+r2t+2b2a1+r2s+2r1t+2ar2+2a1s+2a1t-2r21+s+t, with s=1/α2,t=1/β2, r1≔a1-a2, r2≔a1-a3.Definep1≔cos⁡ϕ,p2≔sin⁡ϕ,q1≔cos⁡θ,q2≔sin⁡θ, and (31)L1≔-6α2p12q22-6β2p22+2p12q12-8α2+β2p1p2-8p1q1,L2≔-40β2p1p2-32(1+sβ2p1p2q1-481+α2p12q1q2-8α24β2+5α2p1p2q2-4β2p22-4p12q12-4α2p12q22-48tα2p12q1q2,L3≔-2β2p22-2α2p12q22-2p12q12,L4≔b2+2b3+3b43p12q12+α2p12q22+β2p22.Then, we have the following proposition which reveals the error between the numerical wavenumber kN and the exact wavenumber k.Proposition 2.
For the 27-point finite difference scheme (7), the numerical wavenumber kN and the exact wavenumber k satisfy(32)kN=k+12a1L1+a2L2+a3L3+L4k3h2+Ok4h3,kh→0.Proof.
Letτ≔kh, and rewrite a,b,c as(33)aτ=cosτcos⁡ϕcos⁡θ,b=cosατcos⁡ϕsin⁡θ,c=cosβτsin⁡ϕ.Hence, M,N can also be considered as two functions with respect to τ. Then, we compute the Taylor expansions of M(τ),1/N(τ) around 0, and after some straightforward manipulations, we obtain(34)Mτ=τ2+τ4a1L1+a2L2+a3L3+Oτ5,τ→0(35)1Nτ=1+τ2L4+Oτ3.τ→0Inserting (34), (35) into (29) yields(36)kNk=1+τ2a1L1+a2L2+a3L3+L4+Oτ3,τ→0Finally, computing again the Taylor expansion of the right-hand side of (36) and replacing τ with kh, we obtain formula (32), which completes the proof.As can be seen from Proposition2, kN converges to k in a second order, which is independent of parameters α and β. This is important since it means that the convergence is robust even if the step sizes in three dimensions are completely different. However, for the 3D rotated difference scheme, it holds(37)kN=μα,βk+νk3h2+Ok4h3,kh→0.where μ(α,β)=1 only for α=β=1 (equidistant sampling). When α≠1 or β≠1 (nonequidistant sampling), μ(α,β)≠1, which means kN does not converge to k, and the rotated scheme fails. Hence, Proposition 2 shows that the new scheme is robust, namely, work effectively for both the situations of equidistant and nonequidistant sampling in different dimensions.


## 3.2. Determination of the Weights

In this subsection, we determine the weights for the new scheme based on minimizing the dispersion error, which is transformed to an optimization problem. We use the least-square method to solve the optimization problem, obtaining a group of optimal weights.For formula (29) in Section 3.1, it is also known as the dispersion relation formula, which is equivalent to the normalized phase velocity [27, 28]. Ideally, kN/k equals 1, which means that the difference scheme is exact. However, for the 2D and 3D situations, it is impossible to achieve this due to the “pollution effect”. Hence, we expect kN/k could be close to 1. For this end, set(38)Ja1,a2,a3,b1,b2,b3,b4;G,ϕ,θ≔kNk-1.When determining the weight parameters a1,a2,a3,b1,b2,b3,b4, (38) is expected to be close to 0. This can be transformed to the following optimization problem:(39)minJa1,a2,a3,b1,b2,b3,b4;G,ϕ,θ∞,s.t.3a1+2a2+a3=1,∑j=14bj=1,Gmin≤G≤Gmax,0≤ϕ≤π2,0≤θ≤π2,where Gmax is usually set to be 400 and Gmin is not smaller than 2 according to the Nyquist sampling limit. To handle (39), let J(a1,a2,a3,b1,b2,b3,b4;G,ϕ,θ)=0, and then obtain(40)A1a1+A2a2+B1b1+B2b2+B3b3=S,where (41)A1≔4bcs+t+4ac1+t+4ab1+s-4a1+r-4bs+r-4ct+r+8r,A2≔8abcr-8abt-4ac1+2s-4bc1+r-4a-4bcs-t-4ct-s-1+4r,B1≔4π2G21-abc,B2≔4π2G213a+b+c-abc,B3≔4π2G213ab+bc+ac-abc,S≔2r-2a+bs+ct-abc4π2G2,with s=1/α2,t=1/β2, r≔1+s+t. Next, we take equidistant samplings for variables G,ϕ and θ in the corresponding intervals Gmin,Gmax,[0,π/2] and [0,π/2] respectively. Assume the numbers of samples for G,ϕ and θ are m,n, and l; then, inserting the samples into (40) yields an overdetermined linear system with unknowns being a1,a2,b1,b2,b3 and size being m×n×l. We use the least-square method to solve the overdetermined linear system and finally obtain a group of optimal weights. For different α and β, the weights vary accordingly. It means that the new difference scheme has a unified form, but equipped with flexible weights, which adaptively vary with the step sizes in three dimensions. In Table 1, we present several groups of weights corresponding to different α and β.Table 1
Weights corresponding to differentα and β.
( α ,  β ) a 1 a 2 a 3 b 1 b 2 b 3 b 4 (1.0, 1.0) 0.1009 0.0011 0.6952 0.3714 0.8111 -0.2928 0.1102 (1.0, 1.5) 0.1154 0.0112 0.6313 0.2614 0.9925 -0.4020 0.1481 (1.5, 1.0) 0.1001 0.0007 0.6984 0.4242 0.6618 -0.1515 0.0656 (1.5, 1.5) 0.1194 0.0061 0.6297 0.3340 0.7907 -0.2159 0.0912 (1.0, 2.0) 0.1522 0.0247 0.4939 0.0274 1.3965 -0.6477 0.2238 (2.0, 1.0) 0.1029 -0.0016 0.6945 0.4275 0.6420 -0.1202 0.0508 (2.0, 2.0) 0.1619 0.0051 0.5040 0.2015 0.9602 -0.2963 0.1346The penalty function in optimization problem (39) is related to the difference between the exact and numerical wavenumbers. Since penalizing all wavenumbers equally, this may not be the most appropriate one. In practice, the more appropriate would be one that penalizes more strongly the error at large wavenumbers, since the cumulative dispersion error is more critical, and the spectral pollution problem is more acute for higher frequencies/wavenumbers. To remedy this, we adopt a flexible strategy for determination of the weights; that is, we estimate the interval IG≔[Gmin,Gmax] by using a priori information. In practice, for example, if the frequency f∈[fmin,fmax] and the velocity v∈[vmin,vmax], then, for a given step size h, we have Gmin≔vmin/hfmax and Gmax≔vmax/hfmin. Hence, the interval IG is flexible according to different cases, and a group of appropriate weights can be obtained for the difference scheme.


## 4. Numerical Experiments

### 4.1. Normalized Phase Velocity Curves

In this subsection, to validate the efficiency and robustness of the new difference scheme, we plot its normalized phase velocity curves associated withkN/k. For comparison, we also present the normalized phase velocity curves for the 27-point rotated difference scheme.According to (29) in Section 3.1, kN/k can be seemed as the function with respect to 1/G if the parameters s a1,a2,a3,b1,b2,b3,b4,α,β,ϕ,θ are given. Plotting kN/k with respect to 1/G obtains the normalized phase velocity curve, which measures the dispersion error. A good normalized phase velocity curve should be close to 1, which indicates a less dispersion error of the corresponding difference scheme. However, the normalized phase velocity curve always deviates gradually from 1 with the decrease of G. From the physical point of view, this is obvious, since the less grid points per wavelength are used, the lower accuracy is obtained.In Figures2 and 3, we present the normalized phase velocity curves for the new 27-point difference scheme and the rotated 27-point difference scheme. For each picture, six curves are plotted, which correspond to the propagation angles (ϕ,θ)=(0,π/4), (π/6,π/3), (π/4,3π/4), (π/3,0), (3π/4,π/6), and (π/4,π/2). Specifically, kN/k is plotted on the vertical y-axis against 1/G on the horizontal x-axis. To validate the robustness of the new difference scheme, we test the cases of (α,β)=(1,1), (1.5,1), (1,1.5), (1,2), (2,1), and (2,1.5), respectively. Thereinto, the case of (α,β)=(1,1) corresponds to the equidistant sampling in all the three dimensions (hx=hy=hz), while other cases correspond to nonequidistant samplings.Figure 2
Normalized phase velocity curves for the new scheme and the rotated scheme with(α,β)=(1,1),(1,1.5),(1.5,1).
(a)
New scheme withα=1,β=1 (b)
Rotated scheme withα=1,β=1 (c)
New scheme withα=1,β=1.5 (d)
Rotated scheme withα=1,β=1.5 (e)
New scheme withα=1.5,β=1 (f)
Rotated scheme withα=1.5,β=1Figure 3
Normalized phase velocity curves for the new scheme and the rotated scheme with(α,β)=(1,2),(2,1),(2,1.5).
(a)
New scheme withα=1,β=2 (b)
Rotated scheme withα=1,β=2 (c)
New scheme withα=2,β=1 (d)
Rotated scheme withα=2,β=1 (e)
New scheme withα=2,β=1.5 (f)
Rotated scheme withα=2,β=1.5As is seen from Figures2 and 3, both the new scheme and the rotated scheme perform efficiently and produce similar dispersion errors when α=β=1. For other situations, the new difference scheme still performs robustly in reducing the dispersion error, which coincides with the Proposition 2. However, the rotated difference scheme fails in this situation. Hence, the new scheme has a decided advantage over the rotated scheme in dealing with the Helmholtz problems of nonequidistant samplings, which is often used in practical applications.We next present some qualitative results. For this end, we compare the numerical dispersion of the proposed 27-point scheme with that of a standard second-order 7-point finite-difference scheme. In Table2, we list the numerical dispersion error kN/k-1 for both schemes corresponding to different parameters α,β,ϕ,θ. Here, we choose G=4, that is, each wavelength possesses 4 grid points in the horizontal direction. As is observed from Table 2, the new 27-point scheme can reduce significantly the numerical dispersion in comparison with the standard 7-point scheme.Table 2
Numerical dispersion of the new scheme versus that of a standard second-order scheme.
Scheme ( α , β ) (0°, 45°) (30°, 60°) (45°, 75°) (75°, 30°) (45°, 90°) 7-point (1.0,1.0) 0.0506 0.0418 0.0474 0.0870 0.0506 27-point (1.0,1.0) 0.0002 0.0006 0.0001 0.0016 0.0009  7-point (2.0,1.0) 0.1191 0.1261 0.1074 0.0871 0.1191 27-point (2.0,1.0) 0.0123 0.0277 0.0009 0.0014 0.0146  7-point (1.0,2.0) 0.0506 0.0601 0.1156 0.3140 0.1191 27-point (1.0,2.0) 0.0118 0.0116 0.0016 0.0355 0.0061  7-point (2.0,2.0) 0.1191 0.1462 0.1806 0.3142 0.1924 27-point (2.0,2.0) 0.0121 0.0364 0.0035 0.0358 0.0637


### 4.2. Numerical Verification of the Convergence Rate

In this subsection, numerical examples are given to verify the convergence order of the new scheme, that is, second order in accuracy.DenoteΔu≔∂2u/∂x2+∂2u/∂y2+∂2u/∂z2. Consider the following 3D Helmholtz equations [29, 30]Example 3.
(42) - Δ u - k 2 u = 3 π 2 - k 2 cos ⁡ π x sin ⁡ π y sin ⁡ π z ,  0 ≤ x , y , z ≤ 1 ,with the Neumann and Dirichlet boundary conditions: ux(0,y,z)=0, ux(1,y,z)=-sin⁡πysin⁡(πz), u(x,0,z)=0, u(x,1,z)=0, u(x,y,0)=0, u(x,y,1)=0. The exact solution is ux,y,z≔cos⁡πxsin⁡πysin⁡πx. In computation, k is chosen to be 10.Example 4.
(43) - Δ u - k 2 u = 3 n 2 - 1 sin ⁡ n k x sin ⁡ n k y sin ⁡ n k z ,  0 ≤ x , y , z ≤ 1 ,with the Dirichlet boundary conditions on all sides of a unit square. When k=mπ,m,n=1,2,⋯, the exact solution is ux,y,z≔sin⁡nkxsin⁡nkysin⁡nkz/k2. In computation, k,m are chosen to be 5π and 1, respectively.We solve the 3D Helmholtz equations (42) and (43) with the new difference scheme and then present the error in Figure 4, where x-axis represents 1/h and y-axis represents the numerical error that is measured by the norm ·∞. It is observed from Figure 4 that the new difference scheme is second order in accuracy. Furthermore, the convergence order is independence of α and β, which further confirms the robustness of the new difference scheme.Figure 4
Numerical convergence order of the new scheme.
(a)
Example3 (b)
Example4


### 4.3. Numerical Examples Related to Geophysical Applications

Finally in this subsection, we test the new difference scheme with solving the 3D Helmholtz equation related to the geophysical applications.We consider the 3D Helmholtz equation in a square domainΛ≔[0,1]×[0,1]×[0,1] with homogeneous medium. To simulate Helmholtz equation in the infinite domain Λ, an artificial boundary condition should be employed. Here, we use the perfectly matched layer (PML, cf. [31, 32]), which has an excellent absorbing performance and produces almost no reflection at the interface. A point source is placed at different locations of the domain Λ, and it generates a Ricker pulse. The simulation results are exhibited in Figure 5, where the real parts of the numerical solutions are visualized with the dimensionless wavenumbers k=30. As can be observed from these figures, there are no obvious dispersion and absorptions, which demonstrate the effectiveness of the new scheme.Figure 5
Visualization of the numerical solutions for the 3D Helmholtz equation.
(a) (b) (c) (d) (e) (f)


## 4.1. Normalized Phase Velocity Curves

In this subsection, to validate the efficiency and robustness of the new difference scheme, we plot its normalized phase velocity curves associated withkN/k. For comparison, we also present the normalized phase velocity curves for the 27-point rotated difference scheme.According to (29) in Section 3.1, kN/k can be seemed as the function with respect to 1/G if the parameters s a1,a2,a3,b1,b2,b3,b4,α,β,ϕ,θ are given. Plotting kN/k with respect to 1/G obtains the normalized phase velocity curve, which measures the dispersion error. A good normalized phase velocity curve should be close to 1, which indicates a less dispersion error of the corresponding difference scheme. However, the normalized phase velocity curve always deviates gradually from 1 with the decrease of G. From the physical point of view, this is obvious, since the less grid points per wavelength are used, the lower accuracy is obtained.In Figures2 and 3, we present the normalized phase velocity curves for the new 27-point difference scheme and the rotated 27-point difference scheme. For each picture, six curves are plotted, which correspond to the propagation angles (ϕ,θ)=(0,π/4), (π/6,π/3), (π/4,3π/4), (π/3,0), (3π/4,π/6), and (π/4,π/2). Specifically, kN/k is plotted on the vertical y-axis against 1/G on the horizontal x-axis. To validate the robustness of the new difference scheme, we test the cases of (α,β)=(1,1), (1.5,1), (1,1.5), (1,2), (2,1), and (2,1.5), respectively. Thereinto, the case of (α,β)=(1,1) corresponds to the equidistant sampling in all the three dimensions (hx=hy=hz), while other cases correspond to nonequidistant samplings.Figure 2
Normalized phase velocity curves for the new scheme and the rotated scheme with(α,β)=(1,1),(1,1.5),(1.5,1).
(a)
New scheme withα=1,β=1 (b)
Rotated scheme withα=1,β=1 (c)
New scheme withα=1,β=1.5 (d)
Rotated scheme withα=1,β=1.5 (e)
New scheme withα=1.5,β=1 (f)
Rotated scheme withα=1.5,β=1Figure 3
Normalized phase velocity curves for the new scheme and the rotated scheme with(α,β)=(1,2),(2,1),(2,1.5).
(a)
New scheme withα=1,β=2 (b)
Rotated scheme withα=1,β=2 (c)
New scheme withα=2,β=1 (d)
Rotated scheme withα=2,β=1 (e)
New scheme withα=2,β=1.5 (f)
Rotated scheme withα=2,β=1.5As is seen from Figures2 and 3, both the new scheme and the rotated scheme perform efficiently and produce similar dispersion errors when α=β=1. For other situations, the new difference scheme still performs robustly in reducing the dispersion error, which coincides with the Proposition 2. However, the rotated difference scheme fails in this situation. Hence, the new scheme has a decided advantage over the rotated scheme in dealing with the Helmholtz problems of nonequidistant samplings, which is often used in practical applications.We next present some qualitative results. For this end, we compare the numerical dispersion of the proposed 27-point scheme with that of a standard second-order 7-point finite-difference scheme. In Table2, we list the numerical dispersion error kN/k-1 for both schemes corresponding to different parameters α,β,ϕ,θ. Here, we choose G=4, that is, each wavelength possesses 4 grid points in the horizontal direction. As is observed from Table 2, the new 27-point scheme can reduce significantly the numerical dispersion in comparison with the standard 7-point scheme.Table 2
Numerical dispersion of the new scheme versus that of a standard second-order scheme.
Scheme ( α , β ) (0°, 45°) (30°, 60°) (45°, 75°) (75°, 30°) (45°, 90°) 7-point (1.0,1.0) 0.0506 0.0418 0.0474 0.0870 0.0506 27-point (1.0,1.0) 0.0002 0.0006 0.0001 0.0016 0.0009  7-point (2.0,1.0) 0.1191 0.1261 0.1074 0.0871 0.1191 27-point (2.0,1.0) 0.0123 0.0277 0.0009 0.0014 0.0146  7-point (1.0,2.0) 0.0506 0.0601 0.1156 0.3140 0.1191 27-point (1.0,2.0) 0.0118 0.0116 0.0016 0.0355 0.0061  7-point (2.0,2.0) 0.1191 0.1462 0.1806 0.3142 0.1924 27-point (2.0,2.0) 0.0121 0.0364 0.0035 0.0358 0.0637


## 4.2. Numerical Verification of the Convergence Rate

In this subsection, numerical examples are given to verify the convergence order of the new scheme, that is, second order in accuracy.DenoteΔu≔∂2u/∂x2+∂2u/∂y2+∂2u/∂z2. Consider the following 3D Helmholtz equations [29, 30]Example 3.
(42) - Δ u - k 2 u = 3 π 2 - k 2 cos ⁡ π x sin ⁡ π y sin ⁡ π z ,  0 ≤ x , y , z ≤ 1 ,with the Neumann and Dirichlet boundary conditions: ux(0,y,z)=0, ux(1,y,z)=-sin⁡πysin⁡(πz), u(x,0,z)=0, u(x,1,z)=0, u(x,y,0)=0, u(x,y,1)=0. The exact solution is ux,y,z≔cos⁡πxsin⁡πysin⁡πx. In computation, k is chosen to be 10.Example 4.
(43) - Δ u - k 2 u = 3 n 2 - 1 sin ⁡ n k x sin ⁡ n k y sin ⁡ n k z ,  0 ≤ x , y , z ≤ 1 ,with the Dirichlet boundary conditions on all sides of a unit square. When k=mπ,m,n=1,2,⋯, the exact solution is ux,y,z≔sin⁡nkxsin⁡nkysin⁡nkz/k2. In computation, k,m are chosen to be 5π and 1, respectively.We solve the 3D Helmholtz equations (42) and (43) with the new difference scheme and then present the error in Figure 4, where x-axis represents 1/h and y-axis represents the numerical error that is measured by the norm ·∞. It is observed from Figure 4 that the new difference scheme is second order in accuracy. Furthermore, the convergence order is independence of α and β, which further confirms the robustness of the new difference scheme.Figure 4
Numerical convergence order of the new scheme.
(a)
Example3 (b)
Example4


## 4.3. Numerical Examples Related to Geophysical Applications

Finally in this subsection, we test the new difference scheme with solving the 3D Helmholtz equation related to the geophysical applications.We consider the 3D Helmholtz equation in a square domainΛ≔[0,1]×[0,1]×[0,1] with homogeneous medium. To simulate Helmholtz equation in the infinite domain Λ, an artificial boundary condition should be employed. Here, we use the perfectly matched layer (PML, cf. [31, 32]), which has an excellent absorbing performance and produces almost no reflection at the interface. A point source is placed at different locations of the domain Λ, and it generates a Ricker pulse. The simulation results are exhibited in Figure 5, where the real parts of the numerical solutions are visualized with the dimensionless wavenumbers k=30. As can be observed from these figures, there are no obvious dispersion and absorptions, which demonstrate the effectiveness of the new scheme.Figure 5
Visualization of the numerical solutions for the 3D Helmholtz equation.
(a) (b) (c) (d) (e) (f)


## 5. Conclusions

For solving the 3D Helmholtz equation, the traditional 7-point central difference scheme gives a poor performance due to its low accuracy. To improve the accuracy, the rotated difference scheme based on staggered-grids was developed, with optimal coefficients obtained by minimizing the numerical dispersion. However, the rotated scheme is rather complex and cannot handle the problems with nonequidistant sampling. Hence, it deteriorates seriously when the step sizes are not equal in all the three dimensions. To overcome this problem, we propose a new optimal 27-point finite difference scheme in this paper. We achieve this by discretizing the Laplacian with a generalized central difference formula, which employs four groups of neighboring grid points, possessing 27 grid points in total. Each group is involved in the difference scheme by the manner of weighted combination. For the discretization of the zeroth term, we use a weighted average of another four groups of grid points, which also possess a total of 27 points. Finally, the optimal weights are determined by minimizing the dispersion error, which is transformed to solving an optimization problem with the least-square method. The new scheme performs robustly for both the problems of equidistant and nonequidistant sampling. We prove this theoretically as well as numerically. Furthermore, the new scheme is much simpler since it is rotation-free. Numerical examples demonstrate the efficiency of the new scheme.



---

*Source: 8532408-2019-07-22.xml*