Research Article  Open Access
A Robust Optimal Finite Difference Scheme for the ThreeDimensional Helmholtz Equation
Abstract
We propose a robust optimal 27point finite difference scheme for the Helmholtz equation in threedimensional 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 zerothorder 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 leastsquare method. In comparison with the rotated difference scheme based on a staggeredgrid 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.
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 oilgas exploration.
All of the numerical methods suffer from the socalled “pollution effect” [8], which leads to the dispersion (phase) errors. In the cases of highdimensional 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 9point difference scheme, based on the staggeredgrid method [10]. The rotated 9point scheme combined the classical Cartesian coordinate system and its rotated systems to discretize the twodimensional (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 9point scheme obtained a good effect and was extended to other situations (cf. [11–13]). In [11], Shin et al. constructed a 2D 25point rotated difference scheme with determining the optimal weights by the singularvalue decomposition method. In [12], Hustedt et al. developed a group of new finite difference methods based on staggeredgrid stencils of various accuracy and grid rotation strategies. For the threedimensional (3D) Helmholtz Equation, [13] employed staggeredgrid stencils on 7 rotated coordinate systems to approximate the differential operators, which is rather complicated. Motivated by the rotated schemes, many weightedaveraging difference schemes were designed by combining the central secondorder 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 righthand 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 27point difference scheme, which remains weighted and second order in accuracy but is rotationfree. 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 zerothorder term, we employ the weighted average of the 27 points. Both the discretizations of the second derivative operator and the zerothorder term use all of the 27 grid points, which are split into four groups. To determine the weights, the leastsquare 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 finitedifference timedomain (FDTD) simulations. For instance, [21] developed the DRP algorithms to minimize the numerical dispersion error in largescale threedimensional 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 WEFDTD 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 Section 2, we present a robust optimal 27point 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 27Point 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 readswhere denotes the pressure wavefield, is the wavenumber with , are the frequency and velocity, respectively, and denotes a source function.
To describe the new finite difference scheme, we use the network of grid points . Here, , , and , with being the step sizes in the , , and directions, respectively. Note that, for the situation of equidistant sampling in all three directions, we have , while for nonequidistant sampling, not all of them are equal. We next present the 27point finite difference stencil with numbering in Figure 1, where the grid points are numbered as with . Moreover, represents the center point of the stencil, and the others denote the points surrounding it. For simplicity, the discretization of a function at point is denoted by or .
The discretization of the 3D Helmholtz equation comprises two segments, one of which is the discretization of the second derivative operator and the other is the discretization of the zerothorder term . To approximate and , we establish a generalized central difference scheme as follows:where with and parameters being the weights. As is observed, (2)(4) reserve the form of classical central difference scheme; however, instead of using a single gridpoint, 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, is a combination of and , , , , , , , , , , . In the direction, all the 27 points are divided into four groups which correspond to , , , and , respectively. For the  and directions, the situations are similar.
To approximate the term of zero order , we utilize the weighted average of all the 27 points as follows:where parameters are another group of weights. As is seen, all the 27 points are involved in the discretization of , and they are also divided into four groups, which correspond to respectively.
Finally, the discretization of equation (1) at point (0,0,0) is expressed aswhich gives a new 27point finite difference scheme for the 3D Helmholtz equation. Substituting (2), (3), (4), and (6) into (7), the new scheme is rewritten aswhere and are the coefficients with
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 satisfy
Proof. We start with handling the formula (2). Let , , , ; then, we haveLet . Applying the Taylor expansion to each term of (11) and (12) yieldswhere , , and . Substituting (14) and (13) into (2) and applying the Taylor expansion again, we obtainwhereSimilarly, for formulas (3) and (4), we havewhere For the formula (6), denote its second, third, and fourth parts by , and respectively. According to the Taylor expansion, we obtainFinally, substituting (15), (17), (18), (20), (21), and (22) into the lefthand side of the difference scheme (7) leads towithAs is observed from (23), the new difference scheme has a secondorder convergence rate when and .
Assume , . 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 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 wavenumber 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.
Assume , , , which gives flexible step sizes in different dimension. When , it corresponds to equidistant sampling; otherwise, it corresponds to nonequidistant sampling. Moreover, let , , and be a constant wavenumber.
For the 3D Helmholtz equation (1) with , its exact solution iswhich is also known as the classical planewave solution. Here, is the imaginary unit and and are the propagation angles from the axis and axis, respectively. On the difference stencil, the discrete form of (25) is given byLet denote the wavelength; then, , and are the numbers of grid points per wavelength in the directions of , , and axis, respectively. For convenience, choose as a benchmark. It is well known that . Together with , we have . Then, substituting the discrete planewave solution into difference scheme (8) with , we obtainwhereLet denote the numerical wavenumber, and replace with in (). Then, some straightforward manipulations yieldwhere with , , .
Define , and Then, we have the following proposition which reveals the error between the numerical wavenumber and the exact wavenumber .
Proposition 2. For the 27point finite difference scheme (7), the numerical wavenumber and the exact wavenumber satisfy
Proof. Let , and rewrite asHence, can also be considered as two functions with respect to . Then, we compute the Taylor expansions of around , and after some straightforward manipulations, we obtainInserting (34), (35) into (29) yieldsFinally, computing again the Taylor expansion of the righthand side of (36) and replacing with , we obtain formula (32), which completes the proof.
As can be seen from Proposition 2, converges to 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 holdswhere only for (equidistant sampling). When or (nonequidistant sampling), , which means does not converge to , 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 leastsquare 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, equals , 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 could be close to . For this end, setWhen determining the weight parameters , (38) is expected to be close to . This can be transformed to the following optimization problem:where is usually set to be and is not smaller than according to the Nyquist sampling limit. To handle (39), let , and then obtainwhere with , . Next, we take equidistant samplings for variables and in the corresponding intervals and respectively. Assume the numbers of samples for and are , and ; then, inserting the samples into (40) yields an overdetermined linear system with unknowns being and size being . We use the leastsquare 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 .

The 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 by using a priori information. In practice, for example, if the frequency and the velocity , then, for a given step size , we have and . Hence, the interval 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 with . For comparison, we also present the normalized phase velocity curves for the 27point rotated difference scheme.
According to (29) in Section 3.1, can be seemed as the function with respect to if the parameters s are given. Plotting with respect to obtains the normalized phase velocity curve, which measures the dispersion error. A good normalized phase velocity curve should be close to , which indicates a less dispersion error of the corresponding difference scheme. However, the normalized phase velocity curve always deviates gradually from with the decrease of . From the physical point of view, this is obvious, since the less grid points per wavelength are used, the lower accuracy is obtained.
In Figures 2 and 3, we present the normalized phase velocity curves for the new 27point difference scheme and the rotated 27point difference scheme. For each picture, six curves are plotted, which correspond to the propagation angles , , , , , and . Specifically, is plotted on the vertical axis against on the horizontal axis. To validate the robustness of the new difference scheme, we test the cases of , , , , , and , respectively. Thereinto, the case of corresponds to the equidistant sampling in all the three dimensions (), while other cases correspond to nonequidistant samplings.
(a) New scheme with
(b) Rotated scheme with
(c) New scheme with
(d) Rotated scheme with
(e) New scheme with
(f) Rotated scheme with
(a) New scheme with
(b) Rotated scheme with
(c) New scheme with
(d) Rotated scheme with
(e) New scheme with
(f) Rotated scheme with
As is seen from Figures 2 and 3, both the new scheme and the rotated scheme perform efficiently and produce similar dispersion errors when . 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 27point scheme with that of a standard secondorder 7point finitedifference scheme. In Table 2, we list the numerical dispersion error for both schemes corresponding to different parameters . Here, we choose , that is, each wavelength possesses 4 grid points in the horizontal direction. As is observed from Table 2, the new 27point scheme can reduce significantly the numerical dispersion in comparison with the standard 7point scheme.

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 . Consider the following 3D Helmholtz equations [29, 30]
Example 3. with the Neumann and Dirichlet boundary conditions: , , , , , . The exact solution is . In computation, is chosen to be .
Example 4. with the Dirichlet boundary conditions on all sides of a unit square. When , the exact solution is . In computation, are chosen to be and , respectively.
We solve the 3D Helmholtz equations (42) and (43) with the new difference scheme and then present the error in Figure 4, where axis represents and 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.
(a) Example 3
(b) Example 4
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 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 . As can be observed from these figures, there are no obvious dispersion and absorptions, which demonstrate the effectiveness of the new scheme.
(a)
(b)
(c)
(d)
(e)
(f)
5. Conclusions
For solving the 3D Helmholtz equation, the traditional 7point central difference scheme gives a poor performance due to its low accuracy. To improve the accuracy, the rotated difference scheme based on staggeredgrids 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 27point 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 leastsquare 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 rotationfree. Numerical examples demonstrate the efficiency of the new scheme.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study is supported in part by the Natural Science Foundation of China under grant 11701389, by the Characteristic Innovation Project from the Educational Department of Guangdong Province under grant 2018GKTSCX043, by the Research Projects of Shenzhen Institute of Information Technology under grants QN201710 and ZY201714, by the Construct Program of the Key Discipline in Hunan Province, by the opening project of Guangdong Province Key Laboratory of Computational Science at the Sun Yatsen University under grants 2018007, and by Hunan Provincial Natural Science Foundation under grant 2019JJ50395.
References
 D. Gordon, R. Gordon, and E. Turkel, “Compact high order schemes with gradientdirection derivatives for absorbing boundary conditions,” Journal of Computational Physics, vol. 297, pp. 295–315, 2015. View at: Publisher Site  Google Scholar
 T. Wu, Z. Chen, and J. Chen, “Optimal 25point finitedifference subgridding techniques for the 2D helmholtz equation,” Mathematical Problems in Engineering, vol. 2016, Article ID 1719846, 16 pages, 2016. View at: Publisher Site  Google Scholar
 D. Cheng, X. Tan, and T. Zeng, “A dispersion minimizing finite difference scheme for the Helmholtz equation based on pointweighting,” Computers & Mathematics with Applications, vol. 73, no. 11, pp. 2345–2359, 2017. View at: Publisher Site  Google Scholar
 Z. Liu, P. Song, J. Li, J. Li, and X. Zhang, “An optimized implicit finitedifference scheme for the twodimensional Helmholtz equation,” Geophysical Journal International, vol. 202, no. 3, pp. 1805–1826, 2015. View at: Publisher Site  Google Scholar
 I. Babuška, F. Ihlenburg, E. T. Paik, and S. A. Sauter, “A generalized finite element method for solving the Helmholtz equation in two dimensions with minimal pollution,” Computer Methods Applied Mechanics and Engineering, vol. 128, no. 34, pp. 325–359, 1995. View at: Publisher Site  Google Scholar
 H. Chen, P. Lu, and X. Xu, “A robust multilevel method for hybridizable discontinuous Galerkin method for the Helmholtz equation,” Journal of Computational Physics, vol. 264, pp. 133–151, 2014. View at: Publisher Site  Google Scholar
 E. Burman, H. Wu, and L. Zhu, “Linear continuous interior penalty finite element method for Helmholtz equation with high wave number: onedimensional analysis,” Numerical Methods for Partial Differential Equations, vol. 32, no. 5, pp. 1378–1410, 2016. View at: Publisher Site  Google Scholar
 I. M. Babuška and S. A. Sauter, “Is the pollution effect of the FEM avoidable for the helmholtz equation considering high wave numbers?” SIAM Journal on Numerical Analysis, vol. 34, no. 6, pp. 2392–2423, 1997. View at: Publisher Site  Google Scholar
 C.H. Jo, C. Shin, and J. H. Suh, “An optimal 9point, finitedifference, frequencyspace, 2D scalar wave extrapolator,” Geophysics, vol. 61, no. 2, pp. 529–537, 1996. View at: Publisher Site  Google Scholar
 Y. Luo and G. Schuster, “Parsimonious staggered grid finitedifferencing of the wave equation,” Geophysical Research Letters, vol. 17, no. 2, pp. 155–158, 1990. View at: Publisher Site  Google Scholar
 C. Shin and H. Sohn, “A frequencyspace 2D scalar wave extrapolator using extended 25point finitedifference operator,” Geophysics, vol. 63, no. 1, pp. 289–296, 1998. View at: Publisher Site  Google Scholar
 B. Hustedt, S. Operto, and J. Virieux, “Mixedgrid and staggeredgrid finitedifference methods for frequencydomain acoustic wave modelling,” Geophysical Journal International, vol. 157, no. 3, pp. 1269–1296, 2004. View at: Publisher Site  Google Scholar
 S. Operto, J. Virieux, P. Amestoy, J. L’Excellent, L. Giraud, and H. B. Ali, “3D finitedifference frequencydomain modeling of viscoacoustic wave propagation using a massively parallel direct solver: a feasibility study,” Geophysics, vol. 72, no. 5, pp. SM195–SM211, 2007. View at: Publisher Site  Google Scholar
 J. Chen, “A generalized optimal 9point scheme for frequencydomain scalar wave equation,” Journal of Applied Geophysics, vol. 92, pp. 1–7, 2013. View at: Publisher Site  Google Scholar
 D. Cheng, Z. Liu, and T. Wu, “A multigridbased preconditioned solver for the Helmholtz equation with a discretization by 25point difference scheme,” Mathematics and Computers in Simulation, vol. 117, pp. 54–67, 2015. View at: Publisher Site  Google Scholar
 B. GosselinCliche and B. Giroux, “3D frequencydomain finitedifference viscoelasticwave modeling using weighted average 27point operators with optimal coefficients,” Geophysics, vol. 79, no. 3, pp. T169–T188, 2014. View at: Publisher Site  Google Scholar
 E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov, “Compact 2D and 3D Sixth order schemes for the Helmholtz equation with variable wave number,” Journal of Computational Physics, vol. 232, no. 1, pp. 272–287, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 T. Wu, “A dispersion minimizing compact finite difference scheme for the 2D Helmholtz equation,” Journal of Computational and Applied Mathematics, vol. 311, pp. 497–512, 2017. View at: Publisher Site  Google Scholar
 W. L. Miranker, “Difference schemes with best possible truncation error,” Numerische Mathematik, vol. 17, no. 2, pp. 124–142, 1971. View at: Publisher Site  Google Scholar
 C. K. Tam and J. C. Webb, “Dispersionrelationpreserving finite difference schemes for computational acoustics,” Journal of Computational Physics, vol. 107, no. 2, pp. 262–281, 1993. View at: Publisher Site  Google Scholar  MathSciNet
 W. Shumin and F. Teixeira, “Dispersionrelationpreserving FDTD algorithms for largescale threedimensional problems,” IEEE Transactions on Antennas and Propagation, vol. 51, no. 8, pp. 1818–1828, 2003. View at: Publisher Site  Google Scholar
 B. Finkelstein and R. Kastner, “FDTD coeffcient modification schemes of the wave and Maxwells equations for controlling order of accuracy and dispersion errors,” in Proceedings of the IEEE Antennas & Propagation Society International Symposium, IEEE, 2007. View at: Google Scholar
 T. Zygiridis and T. Tsiboukis, “Development of higher order FDTD schemes with controllable dispersion error,” IEEE Transactions on Antennas and Propagation, vol. 53, no. 9, pp. 2952–2960, 2005. View at: Publisher Site  Google Scholar
 S. Wang and F. Teixeira, “Griddispersion error reduction for broadband FDTD electromagnetic simulations,” IEEE Transactions on Magnetics, vol. 40, no. 2, pp. 1440–1443, 2004. View at: Publisher Site  Google Scholar
 B. Finkelstein and R. Kastner, “Finite difference time domain dispersion reduction schemes,” Journal of Computational Physics, vol. 221, no. 1, pp. 422–438, 2007. View at: Publisher Site  Google Scholar
 B. Finkelstein and R. Kastner, “A comprehensive new methodology for formulating FDTD schemes with controlled order of accuracy and dispersion,” IEEE Transactions on Antennas & Propagation, vol. 56, no. 11, pp. 3516–3525, 2008. View at: Publisher Site  Google Scholar  MathSciNet
 L. N. Trefethen, “Group velocity in finite difference schemes,” SIAM Review, vol. 24, no. 2, pp. 113–136, 1982. View at: Publisher Site  Google Scholar
 L. Brillouin, Wave Propagation and Group Velocity, Academic Press, 2013. View at: MathSciNet
 F. Ghaffar, N. Badshah, and S. Islam, “Multigrid method for solution of 3D helmholtz equation based on HOC schemes,” Abstract and Applied Analysis, vol. 2014, Article ID 954658, 14 pages, 2014. View at: Publisher Site  Google Scholar
 G. Sutmann, “Compact finite difference schemes of sixth order for the Helmholtz equation,” Journal of Computational and Applied Mathematics, vol. 203, no. 1, pp. 15–31, 2007. View at: Publisher Site  Google Scholar  MathSciNet
 J. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994. View at: Publisher Site  Google Scholar  MathSciNet
 S. Johnson, Notes on Perfectly Matched Layers (PMLs), Massachusetts Institute of Technology, Massachusetts, USA, 2008.
Copyright
Copyright © 2019 Dongsheng Cheng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.