# An Optimal Fourth-Order Finite Difference Scheme for the Helmholtz Equation Based on the Technique of Matched Interface Boundary

Authors: Dongsheng Cheng; Jianjun Chen; Guangqing Long
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this paper, a new optimal fourth-order 21-point finite difference scheme is proposed to solve the 2D Helmholtz equation numerically, with the technique of matched interface boundary (MIB) utilized to treat boundary problems. For the approximation of Laplacian, two sets of fourth-order difference schemes are derived firstly based on the Taylor formula, with a total of 21 grid points involved. Then, a weighted combination of the two schemes is employed in order to reduce the numerical dispersion, and the weights are determined by minimizing the dispersion. Similarly, for the discretization of the zeroth-order derivative term, a weighted average of all the 21 points is implemented to obtain the fourth-order accuracy. The new scheme is noncompact; hence, it encounters great difficulties in dealing with the boundary conditions, which is crucial to the order of convergence. To tackle this issue, the matched interface boundary (MIB) method is employed and developed, which is originally used to accommodate free edges in the discrete singular convolution analysis. Convergence analysis and dispersion analysis are performed. Numerical examples are given for various boundary conditions, which show that new scheme delivers a fourth order of accuracy and is efficient in reducing the numerical dispersion as well.


---

## Body

## 1. Introduction

The Helmholtz equation, also known as frequency-domain wave equation, is widely applied in many areas of science and engineering, such as in aviation, marine technology, geophysics, and optics. The Helmholtz equation is so important that its numerical simulation has attracted significant research interest. For solving the Helmholtz equation, we mainly have finite element methods (cf. [1–5]) and finite difference methods (see [6–12]). Due to the oscillatory characteristic of waves, all of the numerical methods are affected by the “pollution effect” (see [1, 2]), which cannot be eliminated in two-dimensional (2D) and three-dimensional (3D) situations. With the “pollution effect,” the numerical wavenumber deviates from the original one, which is referred to as the “numerical dispersion” (see [13, 14]). Under such circumstances, more finer grids are required for the numerical methods to ensure the accuracy as the wavenumber increases. In order to suppress the “pollution effect” and reduce the numerical dispersion, many numerical methods have been proposed in the past decades (cf. [1, 2, 15–19]).In recent years, finite difference methods have been frequently used for solving the Helmholtz equation, especially in the engineering fields such as oil-gas exploration. Since the standard 2D 5-point finite difference scheme leads to serious numerical dispersion, many dispersion minimizing finite difference methods were developed. In [9], a rotated 9-point finite difference scheme was proposed by linearly combining the two discretizations of the equation on the classical Cartesian coordinate system and the 45o rotated system and the combination coefficients are derived by minimizing the numerical dispersion. In [10], the rotated 9-point difference scheme was extended to a 25-point scheme. In [6], an optimal 9-point difference scheme is developed, which is consistent with the Helmholtz equation with perfectly matched layer (PML). To improve the robustness, the authors of [7, 16] established robust optimal difference schemes which works efficiently even if the step sizes along different directions are not equal. To reduce the numerical error, higher-order finite difference schemes (cf. [8, 11, 12, 17, 20]) were also commonly used. The higher-order schemes can be grouped into two types; one is compact and other one is noncompact. The compact difference scheme possesses not more than 9 and 27 points for the 2D and 3D cases, respectively. For a compact scheme, it requires the right-hand side term to be smooth enough, which cannot be satisfied in many cases. The noncompact scheme is not as demanding as the compact one, but it encounters difficulties in dealing with the boundary conditions. Hence, convergence order of the noncompact scheme tends to deteriorate when it is not handled properly for the discretization near the boundary.In this text, we propose an optimal noncompact fourth-order 21-point finite difference scheme for solving the 2D Helmholtz equation. To discretize second-order derivative term (Laplacian), we first develop two sets of fourth-order difference schemes based on the Taylor formula, which take up a total of 21 grid points. For the discretization of the zeroth-order derivative term, we also design three sets of fourth-order approximation schemes similarly. Then, by weighting all of the approximation schemes above linearly, we obtain a noncompact fourth-order 21-point finite difference scheme. The weights are determined by minimizing the numerical dispersion, which can be transformed to solving an optimization problem by the least-squares method. In order to prevent the convergence order from deteriorating, we employ the matched interface boundary (MIB, cf. [21–24]) method to handle the boundary condition. The MIB method was originally proposed to accommodate free edges in the discrete singular convolution analysis of beams, plates, and shells. Now, it is often used to solve the problems of curved interfaces for finite difference methods. Here, we extend the MIB technique to treat the discretization near the boundary with Robin and Dirichlet boundary conditions and obtain very satisfactory results. Both theoretical analysis and numerical examples show that the new scheme not only delivers a fourth-order convergence rate but also performs efficiently in reducing the numerical dispersion.The rest of this paper is organized as follows: Section2 is devoted to the construction of an optimal fourth-order 21-point finite difference scheme for the 2D Helmholtz equation and the convergence analysis. Section 3 presents the numerical dispersion analysis and introduces the determination of the weights in the new scheme. Section 4 generalizes the MIB method to treat the discretization near the boundary with Robin and Dirichlet boundary conditions. In Section 5, numerical experiments are given to validate the numerical convergence order and effectiveness of the new scheme. Finally, in Section 6, some conclusions end this paper.


## 2. A Fourth-Order 21-Point Finite Difference Scheme

In this section, we propose a new fourth-order 21-point finite difference scheme for solving the 2D Helmholtz equation and present the convergence analysis.We consider the 2D Helmholtz equation:(1)−Δu−k2u=finΩ,where Δ≔∂2/∂x2+∂2/∂y2 is the Laplacian, k is the wavenumber, u is the unknown, and f is the right-hand side term. For the boundary condition, we usually have the following:Dirichlet boundary condition(2)u=g,in∂Ω.Robin boundary condition(3)iku+∂u∂n=g,in∂Ω.Here, ∂/∂n stands for the outward normal derivative.We next introduce the construction of the new 21-point finite difference scheme. We begin with presenting the 21-point finite difference stencil with numbering in Figure1, where m,n with m,n∈ℕ represents the center point of the stencil, while m+l,n+jl,j=0,±1,±2 and l,j≠0,0 denote the neighboring points of m,n. For convenience, let m,n be identified with mh,nh where h represents the step size. For a known function φ, let φ|m,n≔φmh,nh denote its discretization at point m,n. If φ is unknown, we let φm,n represent its difference approximation at m,n.Figure 1
The 21-point finite difference stencil with numbering.Firstly, we use the classic fourth-order central difference formulas to approximate∂2u/∂x2. For l=0,±1,±2, define(4)φx,hu|m,n+l≔um+1,n+l−2um,n+l+um−1,n+lh2,(5)φx,2hu|m,n≔um+2,n−2um,n+um−2,n4h2,which are the second-order central difference approximations along the vertical direction. With (4) and (5), we continue to define(6)ψ1,xu|m,n≔43φx,hu|m,n−13φx,2hu|m,n,(7)ψ2,xu|m,n≔23φx,hu|m,n+1+φx,hu|m,n−1−16φx,hu|m,n+2+φx,hu|m,n−2+13φx,hu|m,n−φx,2hu|m,n.Then, together with (6) and (7), the approximation of ∂2u/∂x2 is given by(8)ℒxu|m,n≔α1ψ1,xu|m,n+α2ψ2,xu|m,n,where α1 and α2 are the weight parameters and satisfy α1+α2=1. In the same way, the approximations of ∂2u/∂y2 are given by(9)ℒyu|m,n≔α1ψ1,yum,n+α2ψ2,yum,n,where(10)ψ1,yu|m,n≔43φy,hu|m,n−13φy,2hu|m,n.(11)ψ2,yu|m,n≔23φy,hu|m+1,n+φy,hu|m−1,n−16φy,hu|m+2,n+φy,hu|m−2,n+13φy,hu|m,n−φy,2hu|m,n,with(12)φy,hu|m+l,n≔um+l,n+1−2um+l,n+um+l,n−1h2,l=0,±1,±2,φy,2hu|m,n≔um,n+2−2um,n+um,n−24h2.We next discretize the zeroth-order termk2u. For this end, let b1, b2, and b3 be another group of weight parameters which satisfy b1+b2+b3=1. Then, we define(13)ℐ1u|m,n≔b1k2u|m,n,(14)ℐ2,xu|m,n≔b223k2u|m+1,n+k2u|m−1,n−16k2u|m+2,n+k2u|m−2,n,(15)ℐ2,yu|m,n≔b223k2u|m,n+1+k2u|m,n−1−16k2u|m,n+2+k2u|m,n−2,(16)ℐ3,xu|m,n≔b323ℐ3,1,xu|m,n−16ℐ3,2,xu|m,n+13ℐ3,3,xu|m,n,(17)ℐ3,yu|m,n≔b323ℐ3,1,yu|m,n−16ℐ3,2,yu|m,n+13ℐ3,3,yu|m,n,where(18)ℐ3,3,xu|m,n≔12k2u|m+1,n+k2u|m−1,n−k2u|m+2,n+k2u|m−2,n,ℐ3,1,yu|m,n≔12k2u|m+1,n+1+k2u|m+1,n−1+k2u|m−1,n+1+k2u|m−1,n−1,ℐ3,2,yu|m,n≔12k2u|m+2,n+1+k2u|m−2,n+1+k2u|m+2,n−1+k2u|m−2,n−1,ℐ3,3,yu|m,n≔12k2u|m,n+1+k2u|m,n−1−k2u|m,n+2+k2u|m,n−2.With (13)–(17), the approximation of the zeroth-order term k2u is given by(19)ℐu|m,n=ℐ1u|m,n+12ℐ2,xu|m,n+ℐ2,yu|m,n+12ℐ3,xu|m,n+ℐ3,yu|m,n.Finally, we obtain a new 21-point finite difference scheme for 2D Helmholtz (1) by combining (8), (9), and (19) as follows:(20)−ℒxu|m,n−ℒyu|m,n−ℐu|m,n=fm,n.As is seen, two sets of difference schemes are developed to approximate the LaplacianΔ=∂2/∂x2+∂2/∂y2, and they employ 21 grid points. Similarly, there are three sets of discretization schemes for approximating the zeroth-order term k2u, which also use 21 grid points. These schemes are combined linearly to approximate the Helmholtz equation with 5 weight parameters, which can be determined by minimizing the dispersion. Later, we shall present approximation of fourth order in accuracy. A straightforward computation with (6)–(17), the new difference scheme, is further given by(21)T1Um−1,n−2+T2Um,n−2+T1Um+1,n−2+T1Um−2,n−1+T3Um−1,n−1+T4Um,n−1+T3Um+1,n−1+T1Um+2,n−1+T2Um−2,n+T4Um−1,n+T5Um,n+T4Um+1,n+T2Um+2,n+T1Um−2,n+1+T3Um−1,n+1+T4Um,n+1+T3Um+1,n+1+T1Um+2,n+1+T1Um−1,n+2+T2Um,n+2+T1Um+1,n+2=fm,n,where Um+l,n+jl,j=0,±1,±2 and l,j≠±2,±2 denote the unknowns and Tjj=1,2,…,5 are the coefficients with(22)T1=α26h2+b324k2,T2=α1−3α212h2+b212k2,T3=−4α23h2−b33k2,T4=−4α1−3α23h2−b23k2,T5=5α1+α2h2+b1k2.We now present the convergence analysis of the new 21-points difference scheme, which reveals the new scheme delivers a fourth order of accuracy.Proposition 2.
The 21-point finite difference scheme (20) is of fourth order in accuracy.Proof.
We first deal with the formulasψ1,xu|m,n and ψ1,yu|m,n. According to the Taylor formula, we have(23)ψ1,xu|m,n=∂2u∂x2|m,n−h490∂6u∂x6|m,n+Oh6,(24)ψ1,yu|m,n=∂2u∂y2|m,n−h490∂6u∂y6|m,n+Oh6.
Forφx,hu|m,n and φx,2hu|m,n, applying the Taylor formula leads to(25)φx,hu|m,n+j=∂2u∂x2|m,n+j+h212∂4u∂x4|m,n+j+h4360∂6u∂x6|m,n+j+Oh6,(26)φx,2hu|m,n=∂2u∂x2|m,n+h23∂4u∂x4|m,n+2h445∂6u∂x6|m,n+Oh6.
Continue to expand (25) with i=±1,±2 at (m, n), and then substitute (25) and (26) into (7); we obtain(27)ψ2,xum,n=∂2u∂x2|m,n+μh4+Oh6,where(28)μ=−1190∂4u∂x4|m,n+124∂6u∂x6|m,n.
Similarly, forψ2,yu|m,n, we have(29)ψ2,yu|m,n=∂2u∂y2|m,n+υh4+Oh6,with(30)υ=−1190∂4u∂y4|m,n+124∂6u∂y6|m,n.
Thus, substituting (23), (24), (27), and (29) into (8) with α1+α2=1 in mind, we derive(31)−ℒxu|m,n−ℒyu|m,n=−∂2u∂x2−∂2u∂y2+ξh4+Oh6,where(32)ξ≔α2μ+ν−α190∂6u∂x6|m,n+∂6u∂y6|m,n.
We next treat the approximation of zeroth-order term (19). For ℐ2,x and ℐ2,y, it holds(33)ℐ2,xu|m,n=b2k2um,n−h46∂4u∂x4|m,n+Oh6,(34)ℐ2,yu|m,n=b2k2um,n−h46∂4u∂y4|m,n+Oh6.
Forℐ3,xu|m,n, we first expand ℐ3,1,xu|m,n, ℐ3,2,xu|m,n, and ℐ3,3,xu|m,n as follows:(35)ℐ3,1,xu|m,n=k2um,n+1+k2um,n−1+k2h22∂2u∂x2|m,n+1+∂2u∂x2|m,n−1+k2h424∂4u∂x4|m,n+1+∂4u∂x4|m,n−1+Oh6,(36)ℐ3,2,xu|m,n=k2um,n+2+k2um,n−2+k2h22∂2u∂x2|m,n+2+∂2u∂x2|m,n−2+k2h424∂4u∂x4|m,n+2+∂4u∂x4|m,n−2+Oh6,(37)ℐ3,3,xu|m,n=−3k2h22∂2u∂x2|m,n−5k2h48∂4u∂x4|m,n+Oh6.
Substituting (35)–(37) into (16) and applying the Taylor formula again, we obtain(38)ℐ3,xu|m,n=b3k2um,n+μ1h4+Oh6,with(39)μ1=b3k2−724∂4u∂x4m,n+13∂2∂y2∂2∂x2m,n.
Then, by performing the same procedure forℐ3,y, it follows that(40)ℐ3,yu|m,n=b3k2um,n+ν1h4+Oh6,with(41)ν1=b3k2−724∂4u∂y4m,n+13∂2∂x2∂2∂y2m,n.
Thus, substituting (13), (33)–(34), and (38)–(40) into (20) with b1+b2+b3=1 in mind, we obtain(42)ℐu|m,n=um,n+γh4+Oh6,where(43)γ=μ1+ν1−b2k26∂4u∂x4m,n+∂4u∂y4m,n.
Finally, by combining (31) and (42), new 21-point difference scheme (20) is rewritten as(44)−∂2u∂x2|m,n−∂2u∂y2|m,n−k2u|m,n−ξ+γh4+Oh6=fm,n,which shows that the new scheme is of fourth order in accuracy.


## 3. Numerical Dispersion Analysis and Determination of the Weight Parameters

In this section, numerical dispersion is firstly analyzed for the new 21-point difference scheme, which reveals the error relation between the numerical wavenumber and the exact wavenumber. Then, the weight parameters in new difference formula (21) are determined by minimizing the numerical dispersion. The normalized phase velocity curves are also presented to demonstrate the effect of the new scheme in reducing the numerical dispersion.Due to the “pollution effect,” the accuracy of the numerical solution deteriorates with increasing wavenumberk. The “pollution effect” can be expressed as an error between the wavenumber of the numerical solution and the wavenumber of the exact solution (k), which is referred as the numerical dispersion. For a finite difference scheme, only a higher order of convergence is not enough, but a smaller numerical dispersion is also required. For Helmholtz problems in 2D and 3D, the numerical dispersion can never be eliminated but can be reduced. Hence, for 21-point difference scheme (20), we want to minimize the numerical dispersion by implementing a dispersion analysis, via which the weight parameters are determined.We begin with a notationG which denotes the number of grid points per wavelength. According to the relationship of wavenumber, velocity, and frequency of the wave, a straightforward computation leads to kh=2π/G. It is known that the exact solution of the homogeneous Helmholtz equation Δu+k2u=0 is the plane wave ux,y≔e−ikxcosθ+ysinθ, where θ is the propagation angle from the x-axis. To do a dispersion analysis, we substitute the discrete plane wave um,n≔e−ikhmcosθ+nsinθ into new difference formula (21) with fm,n=0 and then obtain the dispersion equation(45)k2b1+b2M0+b3N0=1h2α1M1+N1+α2M2+N2,where(46)M0≔164P1+Q1−P2+Q2,N0≔168P1Q1−P1Q2+P2Q1+P1−P2+Q1−Q2,M1≔−83Q1+Q26+52,N1≔−83P1+P26+52,M2≔23P2Q1−83P1Q1+Q26−23Q1+83P1−23P2+12,N2≔23P1Q2−83P1Q1+P26−23P1+83Q1−23Q2+12,with(47)P1≔coskhsinθ=cos2πGsinθ,P2≔cos2khsinθ=cos4πGsinθ,Q1≔coskhcosθ=cos2πGcosθ,Q2≔cos2khcosθ=cos4πGcosθ.In (45), k is not what it was, it is actually a numerical wavenumber. Hence, we replace it with a new notation kN, and together with α1+α2=1 and b1+b2+b3=1, we have(48)kN=1hM1+M2+α2M2+N2−M1−N11+b2M0−1+b3N0−1.The following proposition reveals the error between the numerical wavenumberkN and the exact wavenumber k.Proposition 3.
For the 21-point finite difference scheme (20), there holds(49)kN=k−ζ1+ζ2k5h4+Ok7h6,kh⟶0.Proof.
Letτ≔kh; then, P1, P2, Q1, and Q2 can be written as functions P1τ=cosτsinθ, Q1τ=cosτcosθ, P2τ=cos2τsinθ, and Q2τ=cos2τcosθ. Similarly, we have functions M0τ,N0τ,M1τ, N1τ,M2τ, and N2τ. Then, with (48), we define(50)f1τ≔M1τ+N1τ+α2M2τ+N2τ−M1τ+N1τ,(51)f2τ≔1+b2M0τ−1+b3N0τ−1.
According to the Taylor expansion, we have(52)P1τ=∑n=03−1nτsinθ2n2n!+Oτ8,(53)P2τ=∑n=03−1n2τsinθ2n2n!+Oτ8,(54)Q1τ=∑n=03−1nτcosθ2n2n!+Oτ8,(55)Q2τ=∑n=03−1n2τcosθ2n2n!+Oτ8.
Substituting (52)–(55) into M0τ, N0τ, M1τ, N1τ,M2τ, and N2τ leads to(56)M0τ=1−τ412sin4θ+cos4θ+Oτ6,(57)N0τ=1−τ46sin4θ+cos4θ+Oτ6,(58)M1τ+N1τ=−83P1+Q1+16P2+Q2+5=τ2−190sin6θ+cos6θτ6+Oτ8,(59)M2τ+N2τ=23P2Q1+P1Q2+2P1+Q1−12P2+Q2−163P1Q1+1=τ2−190sin6θ+cos6θ+16sin4θcos2θ+sin2θcos4θτ6+Oτ8.
Then, substituting (58) and (59) into (50) yields(60)f1τ=α1M1τ+N1τ+α2M2τ+N2τ=τ2−ζ1τ6+Oτ8,with α1+α2=1 and(61)ζ1=190sin6θ+cos6θ+α26sin4θcos2θ+sin2θcos4θ.
Similarly, substituting (56) and (57) into (51) leads to(62)f2τ=b1+b2M0τ+b3N0τ=1−ζ2τ2+Oτ,with b1+b2+b3=1 and(63)ζ2=112b2+2b3sin4θ+cos4θ.
Then, applying Taylor expansion to1/f2τ at τ=0 obtains(64)1f2τ=1−ζ2τ4+Oτ6.
According to (48), (60), and (64), it holds(65)kN=1hf1τf2τ=k1−ζ1+ζ2k4h4+Ok6h6.
Finally, applying the Taylor expansion again to1−ζ1+ζ2k4h4+Ok6h6 at τ=kh=0 leads to the conclusion of this proposition.It is seen from the above proposition thatkN approximates k at a rate of fourth order. However, the term ζ1+ζ2k5h4 are dependent on the wavenumber k, and hence, when k is very large, the accuracy of the approximation deteriorates. This is a reflection of the “pollution effect.”We are now ready to determine the weight parametersα1,α2,b1,b2,andb3 for new 21-point scheme (20) by minimizing the numerical dispersion. Since h=2π/Gk, substituting it into (48) yields(66)kNk=G2πM1+M2+α2M2+N2−M1−N11+b2TM−1+b3TN−1,which is called the dispersion relation formula. According to [6, 9, 10], kN/k is considered as the normalized phase velocity (cf. [25]), which is used to measure the numerical dispersion. Ideally, it is expected that kN/k would be 1. However, this is impossible. As an alternative, we can choose a group of weight parameters α2,b2,andb3 such that kN approximates k as closely as possible, which is equivalent to minimizing the numerical dispersion. For this end, let(67)Jα2,b2,b3;G,θ≔kNk−1=G2πM1+M2+α2M2+N2−M1−N11+b2TM−1+b3TN−1−1.Then, minimizing the numerical dispersion is transformed to solve the optimization problem as follows:(68)α2,b2,b3=argminG,θJα2,b2,b3,G,θ∞,Gmin≤G≤Gmax,0≤θ≤π4.Here, the maximum value ofθ is set to be π/4 due to the symmetry of the plane wave. The values of Gmin and Gmax shall be set according to the wave’s velocity, frequency, and the step size in practice. To solve optimization problem (68), let Jα2,b2,b3,G,θ=0 and derive(69)α2M2+N2−M1+N1−4π2G2b2M0−1+b3N0−1=4π2G2−M1+N1.Then, an overdetermined linear system is obtained from (69) by applying equidistant sampling for Gandθ in intervals Gmin,Gmax and 0,π/4, respectively. Finally, we derive a group of weight parameters α2,b2,andb3 by solving the overdetermined linear system with the least-squares method. In Table 1, different groups of weight parameters are presented for different intervals of G.Table 1
Weight parameters for different intervals of.G
[2, 3][3, 4][4, 5][5, 6][6, 8][8, 10][10, 400]α20.315 410.148 620.108 260.092 160.083 450.075 930.072 44b20.328 710.235 540.209 510.198 470.192 190.186 790.184 15b3−0.02270−0.02616−0.02526−0.02464−0.02418−0.02382−0.02359In Figures2(a) and 2(b), we present the normalized phase velocity curves for the new fourth-order 21-point (21p) difference scheme. As a contrast, the normalized phase velocity curves for the optional second-order 9-point (9p) scheme (cf. [6]) are also given in Figures 2(c) and 2(d). In each picture of Figure 2, we plot four normalized phase velocity curves for different propagation angles which include 0o, 15o, 30o, and 45o. The x-axis denotes 1/G, and the y-axis denotes the normalized phase velocity kN/k. Figures 2(a) and 2(c) correspond to case of G∈4,400, while Figure 2(b) and 2(d) correspond to case of G∈2,400. As can be seen, the new fourth-order 21-point difference scheme is efficient in reducing the numerical dispersion and performs much better than the optional second-order 9-point scheme. We also observe that a smaller G gives rise to a larger numerical dispersion. According to the Nyquist sampling limit, G cannot be less than 2.Figure 2
Normalized phase velocity curves. (a) The 21p scheme forG∈4,400. (b) The 21p scheme for G∈2,400. (c) The 9p scheme for G∈4,400. (d) The 9p scheme for G∈2,400.
(a)(b)(c)(d)


## 4. Matched Interface Boundary Technique for the New Finite Difference Method

In this section, we extend the matched interface boundary (MIB) method to treat the boundary condition, which is of great importance for the new difference scheme to maintain a fourth-order accuracy in practice.We point out that the newly proposed difference scheme here is a noncompact one, since its stencil is wide and possesses 21 grid points. Consequently, it encounters difficulty in dealing with boundary conditions, because the wide difference stencil will refer to grid points outside the computational domain. Thus, special numerical treatments are required near boundaries. For this purpose, we here turn to the MIB method, which is originally used to accommodate free edges in the discrete singular convolution analysis [24] and later generalized to other situations (cf. [21, 23]). Based on the idea of MIB method, we shall establish a set of techniques to treat the boundary problems for the new 21-point difference scheme.


### 4.1. A Review of the Matched Interface Boundary Method

In order to describe clearly the boundary treatmenr of the new difference scheme, we first review briefly the MIB method.For convenience, we consider the 1D Robin boundary condition:(70)iku+∂u∂n=g,where i is the imaginary unit and ∂u/∂n is the outer unit normal derivative. The main idea of the MIB method is to construct fictitious points outside the boundaries such that the wide difference stencil can be used successfully on or near the boundary points. Then, difference stencil should be modified accordingly, that is, to determine the function values on the fictitious points.Assume there areM fictitious points outside the boundary and L+1 grid points inside the computational domain (inner points). Denote the function values on fictitious points by fj with j=1,2,…,M, which are unknown. Denote the function values on inner points by uj with j=0,1,2,…,L, which are known. Here, u0 is a function value on the boundary point, which is represented by x=a. To make it clear, we visualize the fictitious and inner points with their function values in Figure 3(a). Next, we seek for a high-order approach to represent fi in terms of uj by means of discretizing boundary condition (70). For the MIB method, it determines the fictitious values fi iteratively by repeatedly matching the boundary condition across the boundary.Figure 3
Illustration of determining function values on the fictitious point with the MIB method. (a)M fictitious points and L inner points near the boundary, (b) the first step to determine f1, (c) the second step to determine f2, and (d) the remaining steps to determine f3,f4,…,fM.
(a)(b)(c)(d)The first step is to gain the function valuef1 at the fictitious point adjacent to the boundary point x=a (Figure 3(b)). In order to achieve high-order accuracy for treating the boundary condition, one-sided difference approximations (cf. [22]) are employed, which involve L+1 grid points on the inner side of the boundary. Thus, boundary condition (70) is approximated as(71)iku0+C2,11f1+∑j=2L+2C2,j1uj−2=g,where C2,j1 are the one-sided difference weights to approximate the first derivative at x=a by using f1,u0,u1,…,uL and(72)C2,j1=dC2,jxdx|x=x2,C2,jx=∏k=S1,k≠jk=S2x−xkxj−xk,with S1andS2 being the summation limits. Here, S1 is the index of the fictitious point corresponding to f1, while S2 is the index of the inner point corresponding to uL. Assume S1=1, then S2=L+2 and xjj=1,2,…,L+2 are the grid-point sequence on which the values are f1, u0, u1, ⋯, uL, respectively. It is observed that the first subscript of C2,j1 is 2, because the boundary point x=a is the second point in the present stencil. According to (72), f1 is solved since u0,u1,…,uL are known. This means f1 can be represented by a linear combination of u0,u1,…,uL.The second step is to gain the function valuef2. After the first step, f1 is known and then we continue to determine f2 on one more fictitious point as shown in Figure 3(c). Similarly, it is accomplished by discretizing the same boundary condition again, but with one new fictitious point. Boundary condition (70) is now approximated as(73)iku0+C3,11f2+C3,21f1+∑j=3L+3C3,j1uj−3=g,where C3,j1 are the one-sided difference weights to approximate first derivative at u0 by using f2,f1,u0,u1, ⋯,uL. In this case, the first subscript of C3,j1 is 3, since x=a is the third point in the present stencil. Then, f2 can be solved from (73), which is also represented in terms of u0,u1,…,uL. For the remaining function values f3,f4,…,fM, they are determined in the same way by repeatedly matching the boundary condition, as can be seen in Figure 3(d).For the Dirichlet boundary condition, the MIB procedure works similarly, but the boundary condition needs a slight modification. Specifically, we derive a new boundary condition containing derivatives based on the Dirichlet boundary condition and the governing equation. Consider the 1D Helmholtz equation with the Dirichlet boundary condition on the left boundary:(74)−uxx−k2u=−1,x∈0,1,ux=0,x=0.As is seen, atx=0, it holds u0=0 and −uxx0−k2u0=−1. Hence, substituting u0=0 into the governing equation leads to a new boundary condition:(75)uxx0=1,x=0.Then, MIB technique can be applied on formula (75) to determine the function values on the fictitious points.


### 4.2. MIB Technique for the New 21-Point Difference Scheme

In this section, we develop the MIB technique for the new 21-point finite difference scheme.For the 21-point scheme, the application of the MIB method shall be more complicated, because it possesses a wide difference stencil. Note that we consider a rectangular domain in this paper. Hence, the fictitious points are distributed around the rectangular domain. For a typical situation, the fictitious points can be divided into three groups, which are along the horizontal, vertical, and diagonal directions, respectively. Figure4 presents fictitious points of different types, where the solid dots represent inner points, fm−2,n+jj=0,1 and fm−1,n+jj=0,1,2 correspond to the fictitious points in the horizontal area, fm+j,n−1j=0,1 and fm+j,n−1j=0,1,2 correspond to the fictitious points in the vertical area, and fm−1,n−2, fm−2,n−1, and fm−1,n−1 correspond to the fictitious points in the diagonal area.Figure 4
Illustration of the fictitious points for the new 21-point difference scheme along the horizontal, vertical, and diagonal directions, respectively.We next introduce the determination of the function values on different types of fictitious points with the Robin boundary condition. For the Dirichlet boundary condition, the procedure is similar as discussed in the previous section. For a rectangular domain0,1×0,1, consider the typical situation of the Robin boundary condition as follows:(76)iku+∂u∂n=g1onΓ1≔x,y|x=0,y∈0,1,(77)iku+∂u∂n=g2onΓ2≔x,y|y=0,x∈0,1.We now determine the fictitious valuesfi iteratively for the new finite difference scheme by repeatedly matching the Robin boundary conditions across the boundaries Γ1 and Γ2.Firstly, we solve the function valuesfm−1,n+jj=0,1,2 and fm+j,n−1j=0,1,2 by(78)iku0+C2,11fm−1,n+l+∑j=2L+2C2,j1um+j−2,n+l=g1,(79)iku0+C2,11fm+l,n−1+∑j=2L+2C2,j1um+l,n+j−2=g2,with l=0,1,…,L. Figure 5(a) is the illustration of the MIB iterative procedure for this step.Figure 5
Illustration of the MIB iterative procedure for the new 21-point scheme. (a) The first step. (b) The second step.
(a)(b)The second step is to determine the function valuesfm−2,n+jj=0,1 and fm+j,n−2j=0,1 by the same iterative procedures with(80)iku0+C3,11fm−2,n+l′+C3,21fm−1,n+l′+∑j=3L+3C3,j1um+j−3,n+l′=g1,(81)iku0+C3,11fm+l′,n−2+C3,21fm+l′,n−1+∑j=3L+3C3,j1um+l′,n+j−3=g2,where l′=1,2 and the fm−1,n+l′ and fm+l′,n−1 had already been obtained according to (78) and (79). The illustration of the iterative procedures for this step is given in Figure 5(b).Thus, the function values are derived on fictitious points in the horizontal and vertical areas. Next, we determine the function values on fictitious points in the diagonal area. Forfm−1,n−1, it is determined by(82)fm−1,n−1≔12f˜m−1,n−1+f^m−1,n−1,where f˜m−1,n−1 and f^m−1,n−1 are solved from(83)iku0+C2,11f˜m−1,n−1+∑j=2L+2C2,j1fm+j−2,n−1=g1,(84)iku0+C2,11f^m−1,n−1+∑j=2L+2C2,j1fm−1,n+j−2=g2.As is seen,fm−1,n−1 can be obtained through the boundary condition (76) or (77) individually. However, taking into account of robustness, we here compute an average of f˜m−1,n−1 and f^m−1,n−1. Finally, for the remaining function values fm−2,n−1 and fm−1,n−2, they are obtained by(85)iku0+C3,11fm−2,n−1+C3,21fm−1,n−1+∑j=3L+3C3,j1fm+j−3,n−1=g1,(86)iku0+C3,11fm−1,n−2+C3,21fm−1,n−1+∑j=3L+3C3,j1fm−1,n+j−3=g2.The illustration of those two steps is presented in Figure6. Thus, in practice, the 21-point difference scheme can be used successfully near the boundary.Figure 6
Illustration of the MIB iterative procedure for determining the function values on fictitious points in the diagonal area. (a)fm−1,n−1. (b) fm−2,n−1 and fm−1,n−2.
(a)(b)


## 4.1. A Review of the Matched Interface Boundary Method

In order to describe clearly the boundary treatmenr of the new difference scheme, we first review briefly the MIB method.For convenience, we consider the 1D Robin boundary condition:(70)iku+∂u∂n=g,where i is the imaginary unit and ∂u/∂n is the outer unit normal derivative. The main idea of the MIB method is to construct fictitious points outside the boundaries such that the wide difference stencil can be used successfully on or near the boundary points. Then, difference stencil should be modified accordingly, that is, to determine the function values on the fictitious points.Assume there areM fictitious points outside the boundary and L+1 grid points inside the computational domain (inner points). Denote the function values on fictitious points by fj with j=1,2,…,M, which are unknown. Denote the function values on inner points by uj with j=0,1,2,…,L, which are known. Here, u0 is a function value on the boundary point, which is represented by x=a. To make it clear, we visualize the fictitious and inner points with their function values in Figure 3(a). Next, we seek for a high-order approach to represent fi in terms of uj by means of discretizing boundary condition (70). For the MIB method, it determines the fictitious values fi iteratively by repeatedly matching the boundary condition across the boundary.Figure 3
Illustration of determining function values on the fictitious point with the MIB method. (a)M fictitious points and L inner points near the boundary, (b) the first step to determine f1, (c) the second step to determine f2, and (d) the remaining steps to determine f3,f4,…,fM.
(a)(b)(c)(d)The first step is to gain the function valuef1 at the fictitious point adjacent to the boundary point x=a (Figure 3(b)). In order to achieve high-order accuracy for treating the boundary condition, one-sided difference approximations (cf. [22]) are employed, which involve L+1 grid points on the inner side of the boundary. Thus, boundary condition (70) is approximated as(71)iku0+C2,11f1+∑j=2L+2C2,j1uj−2=g,where C2,j1 are the one-sided difference weights to approximate the first derivative at x=a by using f1,u0,u1,…,uL and(72)C2,j1=dC2,jxdx|x=x2,C2,jx=∏k=S1,k≠jk=S2x−xkxj−xk,with S1andS2 being the summation limits. Here, S1 is the index of the fictitious point corresponding to f1, while S2 is the index of the inner point corresponding to uL. Assume S1=1, then S2=L+2 and xjj=1,2,…,L+2 are the grid-point sequence on which the values are f1, u0, u1, ⋯, uL, respectively. It is observed that the first subscript of C2,j1 is 2, because the boundary point x=a is the second point in the present stencil. According to (72), f1 is solved since u0,u1,…,uL are known. This means f1 can be represented by a linear combination of u0,u1,…,uL.The second step is to gain the function valuef2. After the first step, f1 is known and then we continue to determine f2 on one more fictitious point as shown in Figure 3(c). Similarly, it is accomplished by discretizing the same boundary condition again, but with one new fictitious point. Boundary condition (70) is now approximated as(73)iku0+C3,11f2+C3,21f1+∑j=3L+3C3,j1uj−3=g,where C3,j1 are the one-sided difference weights to approximate first derivative at u0 by using f2,f1,u0,u1, ⋯,uL. In this case, the first subscript of C3,j1 is 3, since x=a is the third point in the present stencil. Then, f2 can be solved from (73), which is also represented in terms of u0,u1,…,uL. For the remaining function values f3,f4,…,fM, they are determined in the same way by repeatedly matching the boundary condition, as can be seen in Figure 3(d).For the Dirichlet boundary condition, the MIB procedure works similarly, but the boundary condition needs a slight modification. Specifically, we derive a new boundary condition containing derivatives based on the Dirichlet boundary condition and the governing equation. Consider the 1D Helmholtz equation with the Dirichlet boundary condition on the left boundary:(74)−uxx−k2u=−1,x∈0,1,ux=0,x=0.As is seen, atx=0, it holds u0=0 and −uxx0−k2u0=−1. Hence, substituting u0=0 into the governing equation leads to a new boundary condition:(75)uxx0=1,x=0.Then, MIB technique can be applied on formula (75) to determine the function values on the fictitious points.


## 4.2. MIB Technique for the New 21-Point Difference Scheme

In this section, we develop the MIB technique for the new 21-point finite difference scheme.For the 21-point scheme, the application of the MIB method shall be more complicated, because it possesses a wide difference stencil. Note that we consider a rectangular domain in this paper. Hence, the fictitious points are distributed around the rectangular domain. For a typical situation, the fictitious points can be divided into three groups, which are along the horizontal, vertical, and diagonal directions, respectively. Figure4 presents fictitious points of different types, where the solid dots represent inner points, fm−2,n+jj=0,1 and fm−1,n+jj=0,1,2 correspond to the fictitious points in the horizontal area, fm+j,n−1j=0,1 and fm+j,n−1j=0,1,2 correspond to the fictitious points in the vertical area, and fm−1,n−2, fm−2,n−1, and fm−1,n−1 correspond to the fictitious points in the diagonal area.Figure 4
Illustration of the fictitious points for the new 21-point difference scheme along the horizontal, vertical, and diagonal directions, respectively.We next introduce the determination of the function values on different types of fictitious points with the Robin boundary condition. For the Dirichlet boundary condition, the procedure is similar as discussed in the previous section. For a rectangular domain0,1×0,1, consider the typical situation of the Robin boundary condition as follows:(76)iku+∂u∂n=g1onΓ1≔x,y|x=0,y∈0,1,(77)iku+∂u∂n=g2onΓ2≔x,y|y=0,x∈0,1.We now determine the fictitious valuesfi iteratively for the new finite difference scheme by repeatedly matching the Robin boundary conditions across the boundaries Γ1 and Γ2.Firstly, we solve the function valuesfm−1,n+jj=0,1,2 and fm+j,n−1j=0,1,2 by(78)iku0+C2,11fm−1,n+l+∑j=2L+2C2,j1um+j−2,n+l=g1,(79)iku0+C2,11fm+l,n−1+∑j=2L+2C2,j1um+l,n+j−2=g2,with l=0,1,…,L. Figure 5(a) is the illustration of the MIB iterative procedure for this step.Figure 5
Illustration of the MIB iterative procedure for the new 21-point scheme. (a) The first step. (b) The second step.
(a)(b)The second step is to determine the function valuesfm−2,n+jj=0,1 and fm+j,n−2j=0,1 by the same iterative procedures with(80)iku0+C3,11fm−2,n+l′+C3,21fm−1,n+l′+∑j=3L+3C3,j1um+j−3,n+l′=g1,(81)iku0+C3,11fm+l′,n−2+C3,21fm+l′,n−1+∑j=3L+3C3,j1um+l′,n+j−3=g2,where l′=1,2 and the fm−1,n+l′ and fm+l′,n−1 had already been obtained according to (78) and (79). The illustration of the iterative procedures for this step is given in Figure 5(b).Thus, the function values are derived on fictitious points in the horizontal and vertical areas. Next, we determine the function values on fictitious points in the diagonal area. Forfm−1,n−1, it is determined by(82)fm−1,n−1≔12f˜m−1,n−1+f^m−1,n−1,where f˜m−1,n−1 and f^m−1,n−1 are solved from(83)iku0+C2,11f˜m−1,n−1+∑j=2L+2C2,j1fm+j−2,n−1=g1,(84)iku0+C2,11f^m−1,n−1+∑j=2L+2C2,j1fm−1,n+j−2=g2.As is seen,fm−1,n−1 can be obtained through the boundary condition (76) or (77) individually. However, taking into account of robustness, we here compute an average of f˜m−1,n−1 and f^m−1,n−1. Finally, for the remaining function values fm−2,n−1 and fm−1,n−2, they are obtained by(85)iku0+C3,11fm−2,n−1+C3,21fm−1,n−1+∑j=3L+3C3,j1fm+j−3,n−1=g1,(86)iku0+C3,11fm−1,n−2+C3,21fm−1,n−1+∑j=3L+3C3,j1fm−1,n+j−3=g2.The illustration of those two steps is presented in Figure6. Thus, in practice, the 21-point difference scheme can be used successfully near the boundary.Figure 6
Illustration of the MIB iterative procedure for determining the function values on fictitious points in the diagonal area. (a)fm−1,n−1. (b) fm−2,n−1 and fm−1,n−2.
(a)(b)


## 5. Numerical Experiments

In this section, numerical experiments are presented to demonstrate the accuracy and efficiency of the proposed difference scheme. All the experiments are implemented on a PC with Intel(R) Core(TM) i7-9750H running at 2.6 GHz and with 16.0 GB random access memory.We first present an example with the Dirichlet boundary condition.Example 1.
(87)−Δu−k2u=g,in0,12,where g≔2π2−1k2sinπkxsinπky with wavenumber k being a positive integer. The exact solution is sinπkxsinπky.
Table2 demonstrates the numerical convergence and error of the new scheme for solving (87) with different wavenumbers k. Here, C.O. denotes the numerical convergence order and error represents the absolute error in the norm ⋅∞. As is observed, it delivers a fourth order of accuracy for the new 21-point difference scheme equipped with the developed MIB technique. To make it more clear, in Figure 7, we plot the error in the form of −log2error with respect to log2N. Here, N−1 is the number of grid points in each direction. In Figure 7(a), the error curves are almost the same, since the errors for different k are very close. Figure 7(b) is a local part of Figure 7(a).
We next consider the examples with Robin and mixed boundary conditions.Table 2
Numerical convergence and error of the new 21-point difference scheme for Example1.
kh1/81/161/321/641/1281/25610Error1.0 102e − 046.3 772e − 063.9 963e − 072.4 993e − 081.5 623e − 099.7 631e − 11C.O.3.985 63.996 23.999 13.999 74.000 020Error1.0 795e − 046.8 355e − 064.2 868e − 072.6 816e − 081.6 764e − 091.0 477e − 10C.O.3.980 13.994 83.998 73.999 73.999 940Error1.0 936e − 046.9 286e − 064.3 459e − 072.7 186e − 081.6 995e − 091.0 624e − 10C.O.3.980 43.994 93.998 73.999 73.999 880Error1.0 970e − 046.9 508e − 064.3 599e − 072.7 274e − 081.7 050e − 091.0 657e − 10C.O.3.980 23.994 83.998 73.999 73.999 9Figure 7
(a) Error curves for Example1 with k=10,20,40,80. (b) A local part of (a).
(a)(b)Example 2.
(88)−Δu−k2u=0,in0,12,with the Robin boundary condition as follows:(89)iku+∂u∂n=ik−k1eik2yonΓ1≔x,y|x=0,y∈0,1,iku+∂u∂n=ik−k2eik1xonΓ2≔x,y|y=0,x∈0,1,iku+∂u∂n=ik+k1eik1+k2yonΓ3≔x,y|x=1,y∈0,1,iku+∂u∂n=ik+k2eik1x+k2onΓ4≔x,y|y=1,x∈0,1.
The exact solution of (88) is ux,y=eik1x+k2y with k1,k2=kcosθ,sinθ.Example 3.
(90)−Δu−k2u=g,in0,12,where g=−4+2ikx2+y2+2k2+k3ix2y2+k2+k3ixy+k3i with the mixed boundary condition as follows:(91)∂u∂n=ik1eik2y+1+kiyonΓ1≔x,y|x=0,y∈0,1,iku−∂u∂n=ik+k2eik1x+1+kix−k2onΓ2≔x,y|y=0,x∈0,1,iku−∂u∂n=ik−k1eik1+k2y−4+k2y2−k2onΓ3≔x,y|x=1,y∈0,1,u=eik1x+k2+2+kix2+1+kix+kionΓ4≔x,y|y=1,x∈0,1.
The exact solution of (88) is ux,y=eik1x+k2y+2+kix2y2+1+kixy+ki with k1,k2=kcosθ,sinθ.
Tables3 and 4 give the numerical convergence and error of the new difference scheme for Example 2 and Example 3, respectively. Figures 8(a) and 8(b) show the error curves of the new scheme for Example 2 and Example 3. It continues to confirm that the new 21-point scheme together with the MIB technique contributes to a fourth-order accuracy, even if the boundary conditions are more complicated.
In studies [15, 17], two noncompact 25-point difference schemes were proposed for the Helmholtz equation with the PML (perfectly matched layer) boundary condition. Both the two schemes are constructed based on the technique of derivative-weighting, and the former is of second order in accuracy while the latter is of fourth order. In the PML, the components of the wide difference stencil are set to be zero straightforwardly when they refer to grid points outside the computational domain. This treatment performs well due to the property of attenuation of the PML. However, for the Dirichlet, Robin, and mixed boundary conditions, this treatment fails and leads to a serious deterioration of the convergence order. This is exactly what this paper targets at. To tackle this problem, we propose a new 21-point fourth-order scheme and extend the MIB method to handle the boundary problems. The new difference scheme and the MIB method are integrated closely to discretize the equation near the boundary and achieve a good performance according to the results of the numerical experiments. In addition, the 21-point scheme takes fewer calculations than the 25-point schemes, since it possesses less grid points.Table 3
Numerical convergence and error of the new 21-point difference scheme for Example2.
kh1/201/401/801/1601/320k=25,θ=π/8Error9.3 321e − 027.7 664e − 035.5 943e − 043.7 952e − 052.4 960e − 06C.O.3.586 93.795 23.881 73.926 9k=50,θ=π/8Error1.4 485e − 011.1 474e − 027.8 251e − 045.1 567e − 053.0 021e − 06C.O.3.658 13.874 13.923 63.965 8k=50,θ=π/4Error4.4 847e −013.1 189e −022.2 030e −031.5 601e −041.0 280e −05C.O.3.845 93.823 53.819 73.923 8k=100,θ=π/4Error6.9 915e − 015.0 558e − 023.5 499e − 032.5 206e − 041.611e − 05C.O.3.789 63.832 13.815 93.923 6Table 4
Numerical convergence and error of the new 21-point difference scheme for Example3.
kh1/201/401/801/1601/320k=80,θ=π/8Error4.7 699e − 033.6 285e − 042.5 073e − 051.6 853e − 061.0 766e − 07C.O.3.716 53.855 23.895 33.968 5k=160,θ=π/8Error1.2 777e − 019.4 729e − 036.6 850e − 044.5 204e − 052.8 506e − 06C.O.3.753 63.824 83.886 43.987 1k=80,θ=π/4Error7.4 237e − 015.3 995e − 023.8 345e − 032.6 132e − 041.6 859e − 05C.O.3.781 23.815 73.875 13.954 2k=160,θ=π/4Error9.1 458e + 006.8 867e − 014.9 388e − 023.3 781e − 032.2 130e − 04C.O.3.731 23.801 53.869 83.932 1Figure 8
(a) Error in the C-norm for Example2. (b) Error in the C-norm for Example 3.
(a)(b)


## 6. Conclusions

In this paper, a new optimal fourth-order 21-point finite difference scheme is proposed to solve the 2D Helmholtz equation numerically. We focus on two key points. The first one is to construct a noncompact difference scheme based on minimizing the numerical dispersion, so as to improve the numerical accuracy efficiently. In comparison with the compact 9-point scheme, the noncompact scheme employs more grid points; hence, it is more capable of reducing the numerical dispersion. Furthermore, to gain a high-order accuracy, the noncompact scheme is not so demanding on the smoothness of the right-hand side term. However, the noncompact scheme encounters great difficulties in dealing with the boundary conditions, since it has a wide difference stencil. If the discretization of the boundary condition is not treated properly, the high-order scheme will suffer from degeneration in the convergence rate. To tackle this problem, we employ and extend the MIB method for the new 21-point scheme, which is another key point of this text. For different boundary conditions, we develop a set of MIB techniques to modify the difference stencil on the grid points near the boundary. The MIB method efficiently prevents the convergence rate from deteriorating. Numerical examples with different boundary conditions are presented to demonstrate the efficiency of the new difference scheme with the MIB technique. In addition, convergence analysis and dispersion analysis are also performed.



---

*Source: 2539272-2021-12-22.xml*