Abstract
We construct modified forward, backward, and central finite difference schemes, specifically for the Helmholtz equation, by using the Bloch wave property. All of these modified finite difference approximations provide exact solutions at the nodes of the uniform grid for the second derivative present in the Helmholtz equation and the first derivative in the radiation boundary conditions for wave propagation. The most important feature of the modified schemes is that they work for large as well as low wave numbers, without the common requirement of a very fine mesh size. The superiority of the modified finite difference schemes is illustrated with the help of numerical examples by making a comparison with standard finite difference schemes.
1. Introduction
In the present era, the phenomenon of waves propagating is an unavoidable part of human life as it widely covers their domestic, social, commercial, security, and defence purpose needs. Consequently, the accurate propagation of waves has been and more importantly is one of the prime concerns for scientists, engineers, physicists, and mathematicians. Interestingly, many of the physical problems are governed by the time harmonic form of the wave equation, which is also well known as the Helmholtz equation [1, 2]. A few examples are the propagation of water waves in coastal regions, the scattering of waves from an elastic body, and the highly common propagation of sound waves underwater (SONAR). It is not surprising that for the above mentioned problems the analytical methods mostly fail to provide us with a solution and hence the demand for efficient and reliable numerical methods to solve the Helmholtz equation is obvious. Much effort has been invested in this regard and today we have a range of numerical methods for solving the Helmholtz equation, for instance, the finite volume method [3], finite element method [4, 5], finite difference method [6–10], and spectral element method [11–13] which are commonly used to simulate waves with both time independent and dependent natures.
The only challenge linked with the cost of numerical simulations for the Helmholtz equation (or in general for wave propagation problems) is the rule of thumb [2, 14, 15] which demands at least ten elements per wavelength for the accurate propagation of waves (meaning dispersion and dissipation free propagation). This constraint is a serious worry for the wave propagation community because for a practical application the physical domain can have a length that is thousands of times the wavelength of the wave. More importantly, this constraint becomes a real challenge when we are interested in the propagation of waves with a very high wave number, as a high wave number means (more oscillations) a smaller wavelength and for this case the amount of computer resources such as memory and CPU time required by conventional (or standard) schemes is prohibitively large.
Standard finite difference schemes have extensively been used [6–10] with the intention of having a precise answer to the major question “How can dispersion be eliminated or reduced with less computational expense?” for the simulation of wave propagation. However, a serious impediment of the standard schemes developed so far is that they may produce optimal results for wave propagation problems for low wave numbers or with very high computational cost but completely fail for high wave number propagations. Hence, the motivation is obvious for the development of discretization schemes which are less expensive and highly accurate and efficient for the propagation of waves with very high wave numbers.
In this regard, Nehrbass et al. [8], in , gave a new dimension to the second-order central finite difference (CFD) scheme by redefining the usual second-order CFD approximation with a central node of instead of . The major limitation of this work was that it required prior information about the general solution of the Helmholtz equation and it did not provide exact results for problems with nonreflecting boundary conditions. However, the faster convergence of this scheme was of no comparison to the standard CFD scheme especially when the problem is solved with Dirichlet boundary conditions. Recently, in 2010, the idea of Nehrbass et al. was further taken up by Wong and Li in [10]. They proposed a different approach to solve the Helmholtz equation and their construction did not require the use of the general solution, but instead they used the Helmholtz equation itself recursively. In addition to that they also made the first derivative present in the radiation boundary condition exact. The most important feature of this new scheme was that it works for low as well as for very high wave numbers without the common weakness of fine mesh size even with radiation boundary conditions.
The salient features of both schemes are easy (cheap) implementation which only requires the replacement of the central node; the sparse tridiagonal structure of the matrix is preserved as one has in the case of the standard second-order CFD scheme; both schemes provide nodally exact values on the interior grid points for the Helmholtz equation at any wave number.
In this paper, we present a modified approach to have nodally exact approximations of first- and second-order derivatives for the Helmholtz equation on uniform grids, by using the Bloch wave property. This will give us dispersion free results in one dimension and optimal results in higher dimensions. With the help of this alternate approach, one can make finite difference (central as well as forward and backward) approximations of any order exact.
The organization of the paper is as follows. In Sections 2 and 3, modified schemes are presented for the Helmholtz equation and radiation boundary conditions, respectively. In Section 4, the easy (cheap) implementation of modified schemes is discussed. In Section 6, modified schemes are presented for two-dimensional problems. Sections 5 and 6 are devoted to numerical examples for both one- and two-dimensional problems, respectively.
2. Modified Finite Difference Schemes for the Helmholtz Equation
In order to motivate the ideas, we consider the one-dimensional Helmholtz equation [2] where is the wave number.
2.1. Modified Central Finite Difference (CFD) Schemes
The standard second- and fourth-order central finite difference approximations of the second-order derivative at the node of the uniformly spaced grid are as follows [6, 16]: with being the distance between two adjacent nodes of the grid . With these approximations, for (1), we obtain the following stencils: The above finite difference approximations (3) and (4) of the Helmholtz equation admit nontrivial solutions of the form which satisfy the property which is known as the discrete Bloch wave property [4]. Moreover, is the discrete wave number provided that it satisfies or which we get upon writing the above expressions for as a series in .
It is evident from (7) that second- and fourth-order approximations of the second derivative present in the Helmholtz equation provide second- and fourth-order accurate finite difference schemes for the Helmholtz equation. Moreover, the plus signs in front of the leading order error terms signify that the finite difference approximations will lag throughout the domain compared with the exact solution meaning that the discrete wave number is overestimated as is clear from Table 1, whereas, in the case of finite elements, the discrete wave number is underestimated [4, 5].
For dispersion free propagation, one requires that both the exact and discrete waves propagate with the exact wave number; that is, , which is possible only when the product for both second- and fourth-order schemes as is evident from expressions (7). Therefore, one is interested in discretization schemes (such as finite difference schemes) in which the discrete dispersion relation is independent of both the mesh size and the wave number . With this in mind, we modify the standard second-order central finite difference scheme for the Helmholtz equation (3) such that the new scheme provides the exact solution at the nodes of the grid meaning that at all nodes of the grid. For this, we replace with in the discrete Bloch wave property defined in (5) and obtain Using the above property in (3) and (4) gives or or The scheme obtained in (10) using the Bloch wave property is already presented in texts [8, 10] with alternative formulations. Interestingly, the above schemes (10) and (12) lead back to (3) and (4) if we do Taylor series of the middle node coefficients and and keep only the terms up to second order. Furthermore, inserting a nontrivial solution of the form into (10) and (12), we obtain , which means propagation is dispersion free (or equivalently, these modified schemes provide the exact solution at the nodes of the grid). In order to construct exact central finite difference schemes of all orders, we follow an expression given in the book of Cohen [1] and present the following generalized expression: where , , and is given by with the coefficient of the central node being given by It is also evident from Table 1 that higher-order accurate schemes such as (4) provide better accuracy with a smaller number of elements for the Helmholtz equation. However, higher-order schemes require a greater number of stencil points and consequently: the bandwidth of the resulting matrix increases which is computationally more expensive to invert and the number of fictitious nodes increases for the nodes near to the boundary and on the boundary itself. Let us reconsider (10) which is valid at all nodes only when we have Dirichlet boundary conditions applied at the end nodes and . Otherwise, when the above scheme is applied at the end nodes, it gives us the fictitious nodes given by Consequently, in order to have a unique solution of the resulting system, one requires additional information to kill these fictitious nodes. This is more challenging for the fourth-order scheme which will produce fictitious nodes when applied at the boundary nodes and as well as for the nodes and adjacent to the boundary. Therefore, compact schemes are constructed with the aim of retaining high accuracy with a smaller number of stencil points [17]. These discrepancies can be avoided if one can make standard forward and backward finite difference schemes exact which we do in the following section.
2.2. Modified Forward and Backward Finite Difference Schemes
We now consider the first-order forward and backward finite difference approximations of the second derivative given by [16] Inserting the above approximations in (1) and performing simplifications, we get the following stencil forms for the forward and backward schemes, respectively: Using (8) in (19) and (20), we obtain modified forward and backward finite difference schemes given by Similarly, one can make forward and backward schemes of any order exact. In Table 3, we made only those forward and backward schemes which keep only three terms in the finite difference stencil exact. Concluding, one can use any of the presented schemes to have nodally exact solutions for the Helmholtz equation when the problem is posed along with Dirichlet boundary conditions. In the following section, we present modified schemes when the problem has radiation boundary conditions.
3. Modified Finite Difference Schemes for Radiation Boundary Conditions
We now make the first-order derivative involved in the radiation boundary condition given by exact. The second-order central finite difference approximation of is [16] which on inserting into (23) gives the following stencil: Inserting a nontrivial solution of the form results in which is not dispersion free so in order to make it dispersion free we make use of the Bloch wave property (8) and end up with the following form: which satisfies . Hence, the modified central finite difference scheme (27) provides an exact solution at the last node of the grid. This scheme (27) was also obtained by Wong and Li [10] but with a different formulation. Scheme (27) has one fictitous node for which needs to be avoided. So, we are required to solve (27) with (10) to get Instead of performing this extra step, one can construct forward and backward schemes for radiation boundary conditions and for that we consider first-order forward and backward finite difference schemes with the first derivative given by Now, inserting (29) into (23) gives standard forward and backward finite difference schemes Similarly, inserting (29) into (23) together with (8), we obtain modified forward and backward finite difference schemes for radiation boundary conditions (23), given by In Table 4, we have listed coefficients for standard and modified schemes for radiation boundary conditions (23), for central, forward, and backward schemes.
4. A Note on the Implementation of the Modified Schemes
An interesting feature of the modified schemes is their easy implementation as one does not need to write a brand new code, but instead one just needs to replace the coefficient of the th node in the standard schemes with the modified coefficient; the modified schemes keep the same bandwidth structure as one has in the case of standard schemes and therefore adds no cost to the implementation, but one obtains highly accurate results.
5. Numerical Examples for One-Dimensional Problems
In order to illustrate the superiority of the modified schemes, we solve (1) on . The numerical error is measured using the discrete norm, defined by , , with representing the analytical solution and the computed numerical solution. Moreover, denotes the number of grid points in a uniformly spaced grid.
5.1. Dirichlet Boundary Conditions Applied at Both Ends
First of all, we solve (1) on along with Dirichlet boundary conditions, given by Here, we give full systems for modified schemes (10), (21), and (22) and standard finite difference schemes (3), (19), and (20)
In Table 5, the dispersion error for a fixed mesh of size and a broad range of wave numbers from to is given. It is evident that in the case of standard schemes (3), (19), and (20), the dispersion error depends upon the nondimensional wave number (Table 7). As for , all schemes (3), (19), and (20) provide good results whereas they fail for . However, modified schemes (10), (21), and (22) provide highly accurate results even for and are consistent with the dispersion relation which is independent of both the wave number and the mesh size . Fu [18] made an effort and constructed a scheme in which the dispersion error in the leading order error term was only independent of the wave number.
5.2. Dirichlet Boundary Condition at Left End and Radiation Boundary Condition at Right End: Propagating Wave
This time, we solve (1) on along with Dirichlet and radiation boundary conditions applied at the left and right ends, respectively,
For this problem, we use both standard and modified central finite difference schemes (3) and (10) for the approximation of the second derivative present in the Helmholtz equation and use all modified schemes central, forward, and backward for the first derivative present in (23) presented in Section 3. The full systems for the modified and standard finite difference schemes are (Table 2) And now, for radiation boundary conditions,
Once again the superiority of the modified schemes compared with standard schemes is evident from Table 6 even for the case of radiation boundary conditions.
6. Modified Schemes for the Two-Dimensional Helmholtz Equation on Square Meshes
We consider the two-dimensional Helmholtz equation given by with , where and are user-specified constants and is the incident angle. Now, using the following approximations, (37) takes the form Equation (39) is the standard five-point second-order central finite difference scheme for (37). In order to construct an exact scheme, for (37), we now use the following Bloch wave property for the two-dimensional case given by and (39) takes the following form: or which is the required exact central finite difference scheme for (37). Now, in order to avoid duplication, one can construct exact forward and backward finite difference schemes for (37) as explained in the case of the one-dimensional Helmholtz equation in Section 2.2.
7. Numerical Examples for Two-Dimensional Problems
We reconsider (37) defined on the square domain , which has a plane wave solution of the form given by In order to illustrate the effectiveness of the modified scheme in comparison to the standard finite difference scheme, we solve (43) on a square mesh of mesh length for two types of boundary conditions, with Dirichlet boundary conditions applied at all edges and with Dirichlet boundary conditions applied at edges and radiation boundary conditions applied at edges Once again the numerical error is measured using the discrete norm, defined by , , with representing the analytical solution and the computed numerical solution. Moreover, and denote the number of grid points whereas and denote the number of elements in the and directions, respectively.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publishing of this paper.