Abstract
We propose a robust optimal 27-point finite difference scheme for the Helmholtz equation in three-dimensional domain. In each direction, a special central difference scheme with 27 grid points is developed to approximate the second derivative operator. The 27 grid points are divided into four groups, and each group is involved in the difference scheme by the manner of weighted combination. As for the approximation of the zeroth-order term, we use the weighted average of all the 27 points, which are also divided into four groups. Finally, we obtain the optimal weights by minimizing the numerical dispersion with the least-square method. In comparison with the rotated difference scheme based on a staggered-grid method, the new scheme is simpler, more practical, and much more robust. It works efficiently even if the step sizes along different directions are not equal. However, rotated scheme fails in this situation. We also present the convergence analysis and dispersion analysis. Numerical examples demonstrate the effectiveness of the proposed scheme.
1. Introduction
Many of the physical problems are governed by the wave equation, which is also well known as the Helmholtz equation in frequency domain. The Helmholtz equation finds a wide application in many fields of science, engineering, and industry. To solve the Helmholtz equation numerically, finite difference methods [1–3]) and finite element methods [4–7] are frequently employed. For its easy implementation and less computational complexity, the finite difference method is usually preferred, especially in the engineering field such as oil-gas exploration.
All of the numerical methods suffer from the so-called “pollution effect” [8], which leads to the dispersion (phase) errors. In the cases of high-dimensional problems, the “pollution effect” cannot be eliminated; then, the dispersion errors always exist. Consequently, a good numerical solver for the Helmholtz equation should have a good performance in reducing the dispersion errors. For finite difference methods, the classical central difference scheme gives a poor performance with producing large errors. To improve the performance, Jo et al. [9] proposed the rotated 9-point difference scheme, based on the staggered-grid method [10]. The rotated 9-point scheme combined the classical Cartesian coordinate system and its rotated systems to discretize the two-dimensional (2D) Helmholtz equation; that is, the discretization of the equation is a weighted combination of discretization on the Cartesian system and that on the rotated system. The optimal weights were determined by minimizing the dispersion errors. The rotated 9-point scheme obtained a good effect and was extended to other situations (cf. [11–13]). In [11], Shin et al. constructed a 2-D 25-point rotated difference scheme with determining the optimal weights by the singular-value decomposition method. In [12], Hustedt et al. developed a group of new finite difference methods based on staggered-grid stencils of various accuracy and grid rotation strategies. For the three-dimensional (3-D) Helmholtz Equation, [13] employed staggered-grid stencils on 7 rotated coordinate systems to approximate the differential operators, which is rather complicated. Motivated by the rotated schemes, many weighted-averaging difference schemes were designed by combining the central second-order difference scheme at the central node and its neighboring nodes (cf. [14–16]). High order schemes [1, 17, 18] are also established to improve the accuracy; however, they impose strict requirements on the smoothness of the right-hand side.
The rotated difference schemes are not robust in the sense that they fail when the step sizes are not equal in different directions. Hence, the rotated schemes cannot handle the problems with nonequidistant sampling in different directions. Furthermore, the rotated schemes are complex, especially for the 3D Helmholtz equation, which is frequently used in practical applications. To obtain a robust difference scheme for the 3D Helmholtz equation, in this paper, we propose a robust optimal 27-point difference scheme, which remains weighted and second order in accuracy but is rotation-free. To discretize the second derivative operator of the 3D Helmholtz equation, the new scheme utilizes a generalized central difference formula, which is obtained by replacing each of the 7 grid points with a weighted combination of it and its 11 neighboring grid points. To discretize zeroth-order term, we employ the weighted average of the 27 points. Both the discretizations of the second derivative operator and the zeroth-order term use all of the 27 grid points, which are split into four groups. To determine the weights, the least-square method is used to solve the optimization problem of minimizing the dispersion error. In comparison with the rotated scheme, the new scheme is simple, easy to use. More importantly, it is robust enough to handle the problems with nonequidistant sampling, while the rotated scheme fails in this situation. We prove this by theoretical analysis as well as numerical experiments.
We point out that both the new scheme and rotated scheme are similar to the “dispersion relation preserving” (DRP) schemes that were originally introduced in [19] and continued with [20], since all of these methods are based on reducing the numerical dispersion error related to wave problems. However, the DRP schemes are usually suitable for the finite-difference time-domain (FDTD) simulations. For instance, [21] developed the DRP algorithms to minimize the numerical dispersion error in large-scale three-dimensional FDTD simulations. In [22], the DRP method was further generalized to include specification of both dispersion behavior and the orders of accuracy (OoAs) for both the wave equation (WE) and Maxwell’s equations (MEs). Based on DRP techniques, [23] presented higher order FDTD simulations with controllable dispersion error. In [24], a general procedure was introduced to reduce grid dispersion error in broadband FDTD simulations. In [25], a unified methodology for FDTD simulation was presented based on certain modifications of the characteristic equation that accompanies any given discretized version of the WE. In [26], this development was expanded into a comprehensive new methodology for the systematic generation of WE-FDTD schemes tailored to the spectrum of the excitation. In this paper, the proposed method differs from the DRP schemes since it is formulated to solve the Helmholtz equation only in spatial domain.
The remainder of this paper is organized as follows. In Section 2, we present a robust optimal 27-point finite difference scheme for the 3D Helmholtz Equation and prove it to be second order in accuracy. In Section 3, we perform dispersion error analysis and obtain a group of optimal weights for the new scheme. In Section 4, numerical experiments are present to validate the efficiency and effectiveness of the new scheme. Finally, Section 5 concludes this paper.
2. A Robust Optimal 27-Point Finite Difference Scheme
In this section, we propose a new finite difference scheme for the 3D Helmholtz equation. We also present the convergence analysis, which shows that new scheme is second order in accuracy.
The 3D Helmholtz equation 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 27-point 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 zeroth-order 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 grid-point, a set of grid points are employed. The grid points consist of a point and its 11 neighbors in the manner of weighted combination; for instance, 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 27-point 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 left-hand side of the difference scheme (7) leads towithAs is observed from (23), the new difference scheme has a second-order 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 plane-wave 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 plane-wave 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 27-point 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 right-hand 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 least-square method to solve the optimization problem, obtaining a group of optimal weights.
For formula (29) in Section 3.1, it is also known as the dispersion relation formula, which is equivalent to the normalized phase velocity [27, 28]. Ideally, 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 least-square method to solve the overdetermined linear system and finally obtain a group of optimal weights. For different and , the weights vary accordingly. It means that the new difference scheme has a unified form, but equipped with flexible weights, which adaptively vary with the step sizes in three dimensions. In Table 1, we present several groups of weights corresponding to different and .
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 27-point 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 27-point difference scheme and the rotated 27-point 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 27-point scheme with that of a standard second-order 7-point finite-difference 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 27-point scheme can reduce significantly the numerical dispersion in comparison with the standard 7-point 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 7-point central difference scheme gives a poor performance due to its low accuracy. To improve the accuracy, the rotated difference scheme based on staggered-grids was developed, with optimal coefficients obtained by minimizing the numerical dispersion. However, the rotated scheme is rather complex and cannot handle the problems with nonequidistant sampling. Hence, it deteriorates seriously when the step sizes are not equal in all the three dimensions. To overcome this problem, we propose a new optimal 27-point finite difference scheme in this paper. We achieve this by discretizing the Laplacian with a generalized central difference formula, which employs four groups of neighboring grid points, possessing 27 grid points in total. Each group is involved in the difference scheme by the manner of weighted combination. For the discretization of the zeroth term, we use a weighted average of another four groups of grid points, which also possess a total of 27 points. Finally, the optimal weights are determined by minimizing the dispersion error, which is transformed to solving an optimization problem with the least-square method. The new scheme performs robustly for both the problems of equidistant and nonequidistant sampling. We prove this theoretically as well as numerically. Furthermore, the new scheme is much simpler since it is rotation-free. Numerical examples demonstrate the efficiency of the new scheme.
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 Yat-sen University under grants 2018007, and by Hunan Provincial Natural Science Foundation under grant 2019JJ50395.