#### 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.

#### 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 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: Section 2 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:where is the Laplacian, is the wavenumber, is the unknown, and is the right-hand side term. For the boundary condition, we usually have the following: Dirichlet boundary condition Robin boundary conditionHere, 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 Figure 1, where with represents the center point of the stencil, while and denote the neighboring points of . For convenience, let be identified with where represents the step size. For a known function , let denote its discretization at point . If is unknown, we let represent its difference approximation at .

Firstly, we use the classic fourth-order central difference formulas to approximate . For , definewhich are the second-order central difference approximations along the vertical direction. With (4) and (5), we continue to define

Then, together with (6) and (7), the approximation of is given bywhere and are the weight parameters and satisfy . In the same way, the approximations of are given bywherewith

We next discretize the zeroth-order term . For this end, let , , and be another group of weight parameters which satisfy . Then, we definewhere

With (13)–(17), the approximation of the zeroth-order term is given by

Finally, we obtain a new 21-point finite difference scheme for 2D Helmholtz (1) by combining (8), (9), and (19) as follows:

As is seen, two sets of difference schemes are developed to approximate the Laplacian , and they employ 21 grid points. Similarly, there are three sets of discretization schemes for approximating the zeroth-order term , 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 bywhere and denote the unknowns and are the coefficients with

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 and . According to the Taylor formula, we haveFor and , applying the Taylor formula leads toContinue to expand (25) with at (*m*, *n*), and then substitute (25) and (26) into (7); we obtainwhereSimilarly, for , we havewithThus, substituting (23), (24), (27), and (29) into (8) with in mind, we derivewhereWe next treat the approximation of zeroth-order term (19). For and , it holdsFor , we first expand , , and as follows:Substituting (35)–(37) into (16) and applying the Taylor formula again, we obtainwithThen, by performing the same procedure for , it follows thatwithThus, substituting (13), (33)–(34), and (38)–(40) into (20) with in mind, we obtainwhereFinally, by combining (31) and (42), new 21-point difference scheme (20) is rewritten aswhich 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 wavenumber . The “pollution effect” can be expressed as an error between the wavenumber of the numerical solution and the wavenumber of the exact solution (), 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 notation 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 . It is known that the exact solution of the homogeneous Helmholtz equation is the plane wave , where is the propagation angle from the -axis. To do a dispersion analysis, we substitute the discrete plane wave into new difference formula (21) with and then obtain the dispersion equationwherewith

In (45), is not what it was, it is actually a numerical wavenumber. Hence, we replace it with a new notation , and together with and , we have

The following proposition reveals the error between the numerical wavenumber and the exact wavenumber .

Proposition 3. *For the 21-point finite difference scheme (20), there holds*

*Proof. *Let ; then, , , , and can be written as functions , , , and . Similarly, we have functions , , and . Then, with (48), we defineAccording to the Taylor expansion, we haveSubstituting (52)–(55) into , , , , and leads toThen, substituting (58) and (59) into (50) yieldswith andSimilarly, substituting (56) and (57) into (51) leads towith andThen, applying Taylor expansion to at obtainsAccording to (48), (60), and (64), it holdsFinally, applying the Taylor expansion again to at leads to the conclusion of this proposition.

It is seen from the above proposition that approximates at a rate of fourth order. However, the term are dependent on the wavenumber , and hence, when 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 for new 21-point scheme (20) by minimizing the numerical dispersion. Since , substituting it into (48) yieldswhich is called the dispersion relation formula. According to [6, 9, 10], is considered as the normalized phase velocity (cf. [25]), which is used to measure the numerical dispersion. Ideally, it is expected that would be 1. However, this is impossible. As an alternative, we can choose a group of weight parameters such that approximates as closely as possible, which is equivalent to minimizing the numerical dispersion. For this end, let

Then, minimizing the numerical dispersion is transformed to solve the optimization problem as follows:

Here, the maximum value of is set to be due to the symmetry of the plane wave. The values of and shall be set according to the wave’s velocity, frequency, and the step size in practice. To solve optimization problem (68), let and derive

Then, an overdetermined linear system is obtained from (69) by applying equidistant sampling for in intervals and , respectively. Finally, we derive a group of weight parameters 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 .

In Figures 2(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 , , , and . The -axis denotes , and the -axis denotes the normalized phase velocity . Figures 2(a) and 2(c) correspond to case of , while Figure 2(b) and 2(d) correspond to case of . 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 gives rise to a larger numerical dispersion. According to the Nyquist sampling limit, cannot be less than 2.

**(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:where is the imaginary unit and 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 are fictitious points outside the boundary and grid points inside the computational domain (inner points). Denote the function values on fictitious points by with , which are unknown. Denote the function values on inner points by with , which are known. Here, is a function value on the boundary point, which is represented by . 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 in terms of by means of discretizing boundary condition (70). For the MIB method, it determines the fictitious values iteratively by repeatedly matching the boundary condition across the boundary.

**(a)**

**(b)**

**(c)**

**(d)**

The first step is to gain the function value at the fictitious point adjacent to the boundary point (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 grid points on the inner side of the boundary. Thus, boundary condition (70) is approximated aswhere are the one-sided difference weights to approximate the first derivative at by using andwith being the summation limits. Here, is the index of the fictitious point corresponding to , while is the index of the inner point corresponding to . Assume , then and are the grid-point sequence on which the values are , , , , , respectively. It is observed that the first subscript of is 2, because the boundary point is the second point in the present stencil. According to (72), is solved since are known. This means can be represented by a linear combination of .

The second step is to gain the function value . After the first step, is known and then we continue to determine 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 aswhere are the one-sided difference weights to approximate first derivative at by using , . In this case, the first subscript of is 3, since is the third point in the present stencil. Then, can be solved from (73), which is also represented in terms of . For the remaining function values , 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:

As is seen, at , it holds and . Hence, substituting into the governing equation leads to a new boundary condition:

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. Figure 4 presents fictitious points of different types, where the solid dots represent inner points, and correspond to the fictitious points in the horizontal area, and correspond to the fictitious points in the vertical area, and , , and correspond to the fictitious points in the diagonal area.

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 domain , consider the typical situation of the Robin boundary condition as follows:

We now determine the fictitious values iteratively for the new finite difference scheme by repeatedly matching the Robin boundary conditions across the boundaries and .

Firstly, we solve the function values and bywith . Figure 5(a) is the illustration of the MIB iterative procedure for this step.

**(a)**

**(b)**

The second step is to determine the function values and by the same iterative procedures withwhere and the and 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. For , it is determined bywhere and are solved from

As is seen, can be obtained through the boundary condition (76) or (77) individually. However, taking into account of robustness, we here compute an average of and . Finally, for the remaining function values and , they are obtained by

The illustration of those two steps is presented in Figure 6. Thus, in practice, the 21-point difference scheme can be used successfully near the boundary.

**(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. *where with wavenumber being a positive integer. The exact solution is .

Table 2 demonstrates the numerical convergence and error of the new scheme for solving (87) with different wavenumbers . 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 with respect to . Here, 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 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.

**(a)**

**(b)**

*Example 2. *with the Robin boundary condition as follows:The exact solution of (88) is with .

*Example 3. *where with the mixed boundary condition as follows:The exact solution of (88) is with .

Tables 3 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.

**(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.

#### Data Availability

No data were used to support this study.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was financially supported by the Natural Science Foundation of Guangdong Province in China under grant 2020A1515010566, the Natural Science Foundation of Guangxi Province in China under grant AD20238065, and the key project of Guangxi Provincial Natural Science Foundation of China under grants 2017GXNSFDA198014 and 2018GXNSFDA050014.