Journal of Applied Mathematics

Volume 2014 (2014), Article ID 673106, 9 pages

http://dx.doi.org/10.1155/2014/673106

Research Article

Modified Finite Difference Schemes on Uniform Grids for Simulations of the Helmholtz Equation at Any Wave Number

Hafiz Abdul Wajid,¹ Naseer Ahmed,² Hifza Iqbal,³ and Muhammad Sarmad Arshad⁴

¹Department of Mathematics, COMSATS Institute of Information Technology, Lahore 54000, Pakistan

²Department of Mechanical Engineering, Faculty of Engineering, CECOS University of IT and Emerging Sciences, Peshawar 25000, Pakistan

³Department of Mathematics and Statistics, The University of Lahore, Lahore 54000, Pakistan

⁴Department of Mathematics, University College of Engineering, Sciences and Technology, Lahore Leads University, Lahore 54000, Pakistan

Received 11 April 2014; Accepted 20 June 2014; Published 12 August 2014

Academic Editor: Song Cen

Copyright © 2014 Hafiz Abdul Wajid 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.

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

Table 1: Analysis of the number of elements needed for a dispersion error of .

Second-order scheme
Fourth-order scheme

(, ) approx.
(, ) approx.

10
(, )
(, )
10.00008889

(, )
(, )
100.0000970

(, )
(, )
1000.000099

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,

Table 2: The coefficients of the central nodes for modified and standard schemes, where the coefficient for standard schemes can be obtained from the series representations of the central node of the modified scheme.

Order
Modified CFD scheme
Standard CFD scheme

2

4

6

2
with

Table 3: Coefficients of the th node for modified and standard forward and backward schemes for (1), where the coefficients for standard schemes can be recovered from the real part of the series representations of the modified coefficients.

Modified coefficient
Standard coefficient

Modified coefficient
Standard coefficient

Table 4: Coefficients of the th node for modified and standard central and forward and backward schemes for (23), where the coefficients for standard schemes can be recovered from the real part of the series representations of the modified coefficients.

Modified coefficient
Standard coefficient

Central

Scheme

Forward

Scheme

Backward

Scheme

Table 5: Comparison of the dispersion error for standard central (3) and modified central (10), standard forward (19) and modified forward (21), and backward (20) and modified backward (22) finite difference schemes for fixed with varying wave numbers for the case of Dirichlet boundary conditions (32).

Scheme (3)
Scheme (10)

Scheme (3)
Scheme (10)

Scheme (19)
Scheme (21)

Scheme (19)
Scheme (21)

Scheme (20)
Scheme (22)

Scheme (20)
Scheme (22)

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.

Table 6: Comparison of the dispersion error for standard central (3) and modified central (10), standard forward (19) and modified forward (21), and backward (20) and modified backward (22) finite difference schemes for fixed with varying wave numbers for the case of radiation boundary conditions (34).

Scheme (3)
Scheme (10)

Scheme (3)
Scheme (10)

Scheme (19)
Scheme (21)

Scheme (19)
Scheme (21)

Scheme (20)
Scheme (22)

Scheme (20)
Scheme (22)

Table 7: Comparison of the dispersion error at an angle of for standard and modified finite difference schemes for fixed with varying wave numbers for the case of (45) and (46) boundary conditions.

Dirichlet BCs (45)

Dirichlet and radiation BCs (46)

Standard
Modified

Standard
Modified

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.

References

G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Scientific Computation, Springer, Berlin, Germany, 2002.

F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, vol. 132 of Applied Mathematical Sciences, Springer, New York, NY, USA, 1998.

R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, UK, 2002.

M. Ainsworth and H. A. Wajid, “Explicit discrete dispersion relations for the acoustic wave equation in d-dimensions using finite element, spectral element and optimally blended schemes,” in Computer Methods in Mechanics, vol. 1 of Advanced Structured Materials, pp. 3–17, Springer, Berlin, Germany, 2004.

M. Ainsworth and H. A. Wajid, “Optimally blended spectral-finite element scheme for wave propagation and nonstandard reduced integration,” SIAM Journal on Numerical Analysis, vol. 48, no. 1, pp. 346–371, 2010.

S. K. Lele, “Compact finite difference schemes with spectral-like resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992.

M. Nabavi, M. H. K. Siddiqui, and J. Dargahi, “A new 9-point sixth-order accurate compact finite-difference method for the Helmholtz equation,” Journal of Sound and Vibration, vol. 307, no. 3–5, pp. 972–982, 2007.

J. W. Nehrbass, J. O. Jevtic, and R. Lee, “Reducing the phase error for finite-difference methods without increasing the order,” IEEE Transaction on Antennas and Propagation, vol. 46, no. 8, pp. 1194–1201, 1998.

I. Singer and E. Turkel, “High-order finite difference methods for the Helmholtz equation,” Computer Methods in Applied Mechanics and Engineering, vol. 163, no. 1–4, pp. 343–358, 1998.

Y. S. Wong and G. Li, “Exact finite difference schemes for solving Helmholtz equation at any wavenumber,” International Journal of Numerical Analysis and Modeling B, vol. 2, no. 1, pp. 91–108, 2011.

M. Ainsworth and H. A. Wajid, “Dispersive and dissipative behavior of the spectral element method,” SIAM Journal on Numerical Analysis, vol. 47, no. 5, pp. 3910–3937, 2009.

D. Komatitsch and J. Tromp, “Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models, oceans, rotation and self-gravitation,” Geophysical Journal International, vol. 150, no. 1, pp. 303–318, 2002.

D. Komatitsch and J. Tromp, “Spectral-element simulations of global seismic wave propagation - I. Validation,” Geophysical Journal International, vol. 149, no. 2, pp. 390–412, 2002.

F. Ihlenburg and I. Babuška, “Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM,” Computers & Mathematics with Applications, vol. 30, no. 9, pp. 9–37, 1995.

F. Ihlenburg and I. Babuška, “Finite element solution of the Helmholtz equation with high wave number. II. The h-p version of the FEM,” SIAM Journal on Numerical Analysis, vol. 34, no. 1, pp. 315–358, 1997.

E. H. Twizell, Computational Methods for Partial Differential Equations, Ellis Horwood Series: Mathematics and its Applications, Ellis Horwood, Chichester, UK, 1984.

T. W. H. Sheu, L. W. Hsieh, and C. F. Chen, “Development of a three-point sixth-order Helmholtz scheme,” Journal of Computational Acoustics, vol. 16, no. 3, pp. 343–359, 2008.

Y. Fu, “Compact fourth-order finite difference schemes for Helmholtz equation with high wave numbers,” Journal of Computational Mathematics, vol. 26, no. 1, pp. 98–111, 2008.