Development of a new sixth order accurate compact scheme for two and three dimensional Helmholtz equation
Neelesh Kumar, Ritesh Kumar Dubey
TL;DR
This paper introduces a sixth order compact finite difference scheme for the 2D and 3D Helmholtz equation that remains accurate and robust even at large wave numbers, supported by theoretical analysis and numerical tests.
Contribution
A novel sixth order compact scheme for the Helmholtz equation with wave number independent leading error term, enhancing robustness for high-frequency problems.
Findings
Scheme achieves sixth order accuracy in 2D and 3D
Numerical tests confirm robustness at large wave numbers
Theoretical convergence analysis supports numerical results
Abstract
In this work, a new compact sixth order accurate finite difference scheme for the two and three-dimensional Helmholtz equation is presented. The main significance of the proposed scheme is that its sixth order leading truncation error term does not explicitly depend on the associated wave number. This makes the scheme robust to work for the Helmholtz equation even with large wave numbers. The convergence analysis of the new scheme is given. Numerical results for various benchmark test problems are given to support the theoretical estimates. These numerical results confirm the accuracy and robustness of the proposed scheme.
| 512 | 1.0555 | 0.1623 |
|---|---|---|
| 1024 | 0.7434 | 0.0425 |
| 2048 | 1.0047 | 0.0100 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsAdvanced Numerical Methods in Computational Mathematics · Electromagnetic Simulation and Numerical Methods · Differential Equations and Numerical Methods
Development of a new sixth order accurate compact scheme for two and three dimensional Helmholtz equation
Neelesh Kumar
Ritesh Kumar Dubey
Abstract
In this work, a new compact sixth order accurate finite difference scheme for the two and three-dimensional Helmholtz equation is presented. The main significance of the proposed scheme is that its sixth order leading truncation error term does not explicitly depend on the associated wave number. This makes the scheme robust to work for the Helmholtz equation even with large wave numbers. The convergence analysis of the new scheme is given. Numerical results for various benchmark test problems are given to support the theoretical estimates. These numerical results confirm the accuracy and robustness of the proposed scheme.
Keywords: Finite difference methods, Compact schemes, Convergence,
Helmholtz equations, Wave number.
AMS subject classifications: 65N06, 65N12, 65N15, 35J25, 65Z05.
1 Introduction
Consider the boundary value problem governed by the Helmholtz equation
[TABLE]
with some associated boundary conditions, where is the wave number which is constant. The function and the solution value are assumed to be sufficiently smooth and they have required continuous and bounded derivatives. represents a harmonic source term. In this case, it is also assumed that the source function and its required derivatives are known explicitly.
The Helmholtz equation is a time-harmonic solution of the wave equation. There are many physical phenomena which are described by the Helmholtz equation. Some of these include water wave propagation, acoustic wave scattering, electromagnetic wave propagation. There has been growing interest in constructing compact highly accurate difference schemes for solving partial differential equations [1, 2, 3, 4, 5], particularly to solve the Helmholtz equations [6, 7, 8, 9, 10]. Various numerical techniques for solving the boundary value problems modeled by the Helmholtz equation have been developed using different approaches, to name a few of them as the finite-difference methods [6, 11, 12], the boundary element method [8], the finite-element methods [13, 14] and the spectral-element methods [9, 15].
It is observed that the quality of the numerical solution of the Helmholtz problems depends on the wave number. Thus, the numerical solution for the Helmholtz equation (1) is highly oscillatory in the case of large wave number . Bayliss et al. [16], Babuska and Sauter [14] discussed that for a given accuracy, the number of grid points increases with increasing the wave number. Thus for a given mesh, the numerical errors are developed with the wave number. There are many numerical schemes which can provide accurate numerical results for large wave number. On the other hand, compact high order finite difference techniques for the numerical solution of the Helmholtz equation have been used widely since they provide high accuracy with less computational costs. Nabavi et al. [6], Singer and Turkel [17] discussed the compact sixth order finite difference schemes for the Helmholtz equations in two dimensions. A compact sixth order accurate scheme for the three dimensional Helmholtz equation is presented by Sutmann [18] for constant and by Turkel et al. [19] for variable wave number. Fu [20] developed fourth order accurate finite difference schemes with high wave number for sufficiently small , where is a spatial mesh step size. In this work, a new compact sixth order finite difference scheme for (1) is presented which is also suitable for high wave numbers. The main significance of this work is that the leading truncation error of the proposed scheme does not depend explicitly on the wave number . It is shown that the scheme is uniquely solvable. The convergence analysis for the proposed scheme is also discussed. The proposed scheme is compared with the standard sixth order schemes [6, 11]. The resulting scheme is also compact [21] in the sense that it involves only the patch of cells adjacent to a given node in the mesh. We expect that the proposed scheme will be useful for the Helmholtz equation with large wave numbers.
The paper is organized as follows. In Section 2, a new compact sixth-order scheme for (1) is derived. In Section 4, it is proved that the proposed scheme is uniquely solvable. The error bound is also established for the proposed scheme and the standard sixth order scheme [6]. In Section 5, the numerical results for various benchmark test problems are given to justify the accuracy of the new scheme. At the last, we concluded our work in Section 6.
2 Formulation of a new sixth-order accurate approximation
In the past, Nabavi [6] and others [11, 17] discussed compact sixth order approximations for the Helmholtz equations. The leading truncation error term in a given difference scheme depends on the solution value , the source function and the wave number which shows that the quality of the numerical solution depends significantly on the wave number . In this work, our effort is to make the scheme whose truncation error is less relevance to the wave number .
2.1 Two-Dimensional case
Consider the two-dimensional Helmholtz equation
[TABLE]
with some Dirichlet boundary
[TABLE]
where with its boundary . Assuming to be sufficiently smooth, consider the discretization of rectangular domain on a compact stencil with mesh step size and in - & - directions respectively. The discrete grid points along - & - directions are defined as
For ease of notations, we denote the following notations in further derivation
[TABLE]
The operators used to approximate derivatives with a minimum stencil width are known as compact. The compact scheme amongst the various finite difference replacements has received more attention due to a minimal width of stencils and easy computations. In contrast, high-order difference schemes formulated with non-compact stencils yield a higher bandwidth of the iteration matrix which involves large arithmetic operations.
Using linear combination of minimum number of grid points, the following compact operators corresponding to the second derivatives of at the grid point are given by
[TABLE]
Similarly, the composite operator for the mixed derivatives of is defined as
[TABLE]
Then for sufficiently smooth , Taylor series expansion to (4)-(5) gives
[TABLE]
[TABLE]
Let be the discrete value of the solution value which satisfies (2), then substituting the - approximations of derivatives from (6) into (2), we get
[TABLE]
where and is the leading truncation error, and is given by
[TABLE]
For getting sixth-order accurate schemes, we have to find fourth and second order approximation of respectively fourth and sixth derivatives of solution values.
To achieve this, differentiating (2) twice with respect to and solving for fourth derivatives, we get
[TABLE]
Further differentiating (10) twice with respect to
[TABLE]
Also from (2), we have
[TABLE]
Using (7), we get -approximations of (12)
[TABLE]
Substituting (13) into (7), we will get -approximation of as
[TABLE]
Now combining equations in (10) and using (14), we get -approximation of at the grid as
[TABLE]
Again combining both equations in (11) and using (10) and (12), we have
[TABLE]
From (2), (6) and (7), we get -approximations of at the grid
[TABLE]
Substituting equations (15) and (17) into (9), we have
[TABLE]
where,
[TABLE]
Substituting from (18) into (8), we get the standard sixth-order scheme [6] for (2)
[TABLE]
where the leading truncation error term is given by
[TABLE]
From (20), it is clear that depends upon the solution value , the source function and it also explicitly depends on the wave number .
For simplicity, we omit subscripts in and in further derivation as well.
Now from (2) and (10), we have
[TABLE]
Now, differentiating (11) twice w.r.t. , and combining both, we get
[TABLE]
Using (16) and (21) into (22) , we have
[TABLE]
Substituting (12), (21), and (23) into (20)
[TABLE]
where . Using the - approximations of the second order derivatives from (6)-(7) into (24), we have
[TABLE]
Substituting from (25) into (19), we get a new sixth order compact scheme for the three dimensional Helmholtz equation (2)
[TABLE]
where
[TABLE]
And the leading truncation error term to new sixth order scheme is given by
[TABLE]
From (28), It is clear that the truncation error does not explicitly depend on the wave number . It depends only on the solution value , the source function .
It is noted that the source function taken here is known at each grid point. Therefore right side of equation (26) is fully known. If the source function in equation (2) is unknown, then we need at most fourth order accurate approximation of and which can be obtained by enlarging the cell stencil. Then the resultant system of equation can be solved using some direct or iterative methods [22, 23, 24].
2.2 Three-Dimensional case
We consider the Helmholtz equation in three dimensions
[TABLE]
subject to the Dirichlet boundary condition , where is a sufficiently differentiable function and with its boundary . We consider the discretization of on a compact 27-point cell stencil with uniform mesh step size in each direction of the coordinate axis. The internal grid points along -, - & - directions are and be the discrete value of the solution at this grid.
For ease of notations, let we denote some index sets , and . Also we denote some operators for partial derivatives by
[TABLE]
The following standard finite difference compact operators for the second derivatives of at the grid point are given by
[TABLE]
Thus for sufficiently smooth , Taylor series expansion to (30) gives
[TABLE]
From (30), The composite operators for the mixed derivatives of
[TABLE]
Let satisfies (29) at the grid , then from (31), we have
[TABLE]
where and the local truncation error is given by
[TABLE]
In order to get -approximation of , -approximation of , differentiating (29) twice w.r.t. and solving for fourth derivatives, we get
[TABLE]
Again by differentiating forcing function of (29), we have mixed derivatives as
[TABLE]
Further differentiating (35) twice w.r.t. for sixth derivatives
[TABLE]
Combining the equations in (36), then we have
[TABLE]
Using -approximations from (32) into (38), we get
[TABLE]
For - approximations of , we consider
[TABLE]
where, can be obtained by using the Taylor series expansion to (40)
[TABLE]
Adding equations in (40) and using (41), we have
[TABLE]
Now using -approximation from (39) into (42), we get approximation
[TABLE]
Adding equations in (35) and using (43), we have at the grid
[TABLE]
To find , adding the equations in (37) and using (35) and (38)
[TABLE]
Substituting approximations from (31)-(32) into (45), we have approximations of at the grid
[TABLE]
Finally, substituting -approximation of from (44) and -approximation of from (46) into (34), we get sixth order leading truncation error
[TABLE]
where,
[TABLE]
Substituting from (47) into (33), we get sixth order scheme [11]
[TABLE]
where,
[TABLE]
The coefficients of the scheme and its truncation error term are given by
[TABLE]
For the ease of notations, here and below we omit subscripts . Now from (29), we have
[TABLE]
Now, differentiating (37) twice w.r.t. , and combining them, we get
[TABLE]
Again from (29), we have
[TABLE]
Combining equations given in (52), we have
[TABLE]
Adding equations in (51), we get
[TABLE]
Substituting (38),(45),(50),(53), and (54) into (49), we get
[TABLE]
where, , and .
Using the - approximations from (31)-(32) into (55)
[TABLE]
Finally, we substitute (56) into (48), to get another sixth order compact scheme for the three dimensional Helmholtz equation (29)
[TABLE]
where be the sixth order leading truncation error term to new scheme and is given by
[TABLE]
where, the coefficients of the new scheme are given by
[TABLE]
It is noted that the truncation error does not explicitly depend on the wave number . It depends only on the solution value , the source function .
3 High order accurate approximation of Neumann boundary
Dirichlet boundary conditions specify the solution value at each node of the boundary. Therefore, in case of Dirichlet boundary, the difference schemes can be used for all interior grid points. However, Neumann boundary conditions specify the derivative of the solution value at some part of the boundary. In the later case, the compact schemes can not be used straightforward for all interior grid points since the values at some boundary points are not given.
In this section, a discretization technique is developed for a Neumann boundary condition such that the discretization is consistent with the given sixth-order accurate difference schemes. Here, we consider a sixth-order approximation for and for the two- and three-dimensional Helmholtz equations respectively.
3.1 Two-dimensional case
Without loss of generality, we consider a sixth order accurate discretization for a Neumann condition in two dimensions, where the function has required continuous and bounded derivatives.
Consider the uniform mesh discretization of the domain as . For the aim, we cosider further discretization of the domain by adding a row consisting of ghost points , outside the domain.
The central difference second order approximation of first derivative gives
[TABLE]
Using the Taylor series expansion at , we get
[TABLE]
Now by successive differentiation of the Helmholtz equation (2), we get
[TABLE]
Neumann condition with its tangential derivatives gives
[TABLE]
Substituting the above derivatives at into (59), we subsequently get a sixth-order approximation for the Neumann boundary
[TABLE]
Similar discretizations hold for the Neumann conditions , and in other directions. Using the difference schemes at the boundary point , we can eliminate the values at the ghost points.
3.2 Three-dimensional case
Here, we now consider a Neumann condition for the three-dimensional Helmholtz equation, assuming the function has required continuous and bounded derivatives in and .
As in case of two dimensions, we also introduce a further discretization of the domain having ghost points , outside the boundary of the domain. Then the central difference approximation of derivative gives
[TABLE]
Taylor series expansion at gives
[TABLE]
Successive differentiation of the Helmholtz equation (29) gives
[TABLE]
Using the Neumann boundary condition and its tangential derivatives into (64), we have
[TABLE]
Substituting the Neumann boundary and derivatives and at into (63), we get
[TABLE]
Similar discretizations hold for the Neumann conditions in other directions. Using the difference schemes at the boundary point , we can eliminate the values at the ghost points.
4 Convergence analysis
In this section, the theoretical analysis for the proposed scheme is presented. It is shown that the proposed scheme is uniquely solvable for sufficiently small , where is the wave number and is the mesh step size. The bound on the error norms to the proposed scheme (26) is also established. It is proved that as the mesh step size approaches to zero such that is sufficiently small, then the solution of the proposed difference scheme converges to the solution of the corresponding boundary value problem. We consider here a convergence analysis for the proposed scheme (26) to the Helmholtz equation (2)-(3) in two dimensions. For the three dimensional case, the treatment is similar to the former.
Let be the approximate solution to so that . Consider the difference scheme (26) in operator form as
[TABLE]
where
[TABLE]
with the corresponding truncation error given by (28). The coefficients ’s are given in (27). Now we will show that the proposed difference scheme is uniquely solvable as .
4.1 Solvability of the Difference Scheme
Rewriting the equation (67) in terms of discrete stencil points
[TABLE]
where the operators and are same as in (5) and
[TABLE]
The difference equation (68), in the simple matrix form, can be written as
[TABLE]
where is a tri-block-diagonal matrix with the following elements
[TABLE]
[TABLE]
[TABLE]
Now we claim that the graph of the matrix is strongly connected and hence we have the following result.
Lemma 4.1**.**
Let be a matrix with its elements given by (70)-(72). When is sufficiently small, where is the wave number and is the grid length, then the directed graph of the matrix is strongly connected.
Proof.
Consider the given matrix with its elements given by (70)-(72). For , the distinct nonzero values of the elements of the matrix are given by a00=10/3, a10=-2/3, a20=-1/6. Therefore, for sufficiently small , all the blocks of the matrix contain nonzero elements in its lower, upper and main diagonal. For constructing the directed graph of the matrix , consider distinct points in a plane denoted by . Then corresponding to each nonzero element of the matrix , connect the point to by a directed arc, directed from to . The graph thus obtained is called a directed graph of the matrix [25, 26]. Such a directed graph of a matrix is depicted for in Figure 1 and its adjacency matrix is given as
[TABLE]
where . Hence from the Figure, it is clear that for any ordered pair of points , there exist a directed path from to and hence the directed graph of the matrix is strongly connected [25], which proves the Lemma. ∎
We next establish the conditions on row sum of the matrix in the following result.
Lemma 4.2**.**
Let , where , be a matrix with its elements given by (70)-(72). When is sufficiently small, where is the wave number and is the grid length, then
- (a)
, 2. (b)
, and , for at least one .
Proof.
Consider the given matrix with its elements given by (70)-(72), Clearly when , then a00=10/3, a20=-1/6, a10=-2/3. it is also noted that gives the values of all diagonal element of the matrix and give the values of the off-diagonal elements of the matrix which are clearly non-positive. Thus .
Now let denotes the row sum of the matrix . Then, for ,
[TABLE]
[TABLE]
[TABLE]
From the row sum of the matrix given by (73)-(75), it is clear that all the row sum are nonnegative and at least one is positive, which completes the proof. ∎
Next we will show that the difference scheme has unique solution.
Theorem 4.3**.**
Consider the difference scheme (67) as
[TABLE]
where the operator and are given in (67). Let is sufficiently small, where is the wave number, is the mesh size. Then, there exist a unique solution of the difference scheme (76).
Proof.
Since from the Lemma 4.1, the directed graph of the matrix is strongly connected and hence the matrix is irreducible [25]. Further, from Lemma 4.2, all off-diagonal elements of the matrix are non-positive and also all row sum of the irreducible matrix are nonnegative with at least one positive row sum which implies that the matrix is monotone [27]. Thus the matrix is irreducible and monotone and hence inverse of the matrix exist and [28]. Consequently, the difference scheme (76) has unique solution, which proves the Theorem. ∎
4.2 Error estimate of the Scheme
In this section, the error bound for the proposed difference scheme (26) is established. For this first we will show that the proposed scheme satisfies the discrete maximum principle [22], then we will establish the discrete regularity result [22, 29]. Both the results will be used to prove our main Theorem for getting error bound. Before we present some results, we need to mention some notations for norms which will be used in the later contexts.
Let be the exact solution to the problem (2)-(3) and be the solution to the difference scheme (26), where is the discrete grid on rectangle domain with its boundary . We let define the sup-norms of the solution values
[TABLE]
Now it is claimed that the discrete maximum principle is satisfied by the present difference scheme (26) in the following result.
Lemma 4.4**.**
Let , be a function defined on the grid . If
[TABLE]
* on , where is given in (67) and is an interior of the discrete grid . Then, for sufficiently small , the maximum (minimum) of the solution value on lies on the boundary of the grid . Where, is the wave number and is the grid size.*
Proof.
For the sake of convenience, we will show that the solution value cannot have a local maximum in the interior of the domain.
Let , , from (68) we have
[TABLE]
where, for , . The operator and are same as in (5). Now let us assume that the solution value is a local maximum, which implies that
[TABLE]
Keeping fixed, using (79) into (78) we have
[TABLE]
Now using into (80), equation (80) becomes
[TABLE]
It is also clear that , hence inequality (81) will become
[TABLE]
which further implies that . In this same manner, one can also prove that
[TABLE]
Now keeping fixed, again using assumptions (79), into (78), and proceeding as above, we have
[TABLE]
Hence, cannot have a local maximum in the interior of the domain . Therefore, the maximum of the solution value on must lies on the boundary of the domain. Following the above procedure, One can also prove that the minimum of the solution value on must lies on the boundary of the domain in the case of . The proof is complete. ∎
We next have the following result for the present scheme.
Lemma 4.5**.**
Let be a rectangular grid with its interior and boundary . , be a grid function defined on the grid , with . Then for sufficiently small ,
[TABLE]
where is given in (67) and be the maximum norm defined on the interior of the grid .
Proof.
Consider a grid function defined on the interior of the discrete grid by
[TABLE]
By the definition of maximum norm, it is obvious that
[TABLE]
Then from (86)-(87), we get the inequalities
[TABLE]
We next define , then for , we have
[TABLE]
where . Then, we have
[TABLE]
Now since , Thus from first inequalities in (88), we have
[TABLE]
In the same manner, using second inequalities in (88), we can have
[TABLE]
Thus from Lemma 4.4, the maximum of and the minimum of must occur on the boundary . Then, using on , we get
[TABLE]
[TABLE]
[TABLE]
Proceeding as above for the minimum attained on the boundary, one can prove
[TABLE]
[TABLE]
Since , thus from (96), we get the required result
[TABLE]
which proves the Lemma. ∎
Now we consider error estimate for scheme (26) in the following theorem.
For this, let we denote
[TABLE]
which will be used in the later context.
Theorem 4.6**.**
Let , and be a solution to the problem (2)-(3) and be the solution to the present difference scheme (26). Then for sufficiently small , there exist the following error estimate for the present scheme
[TABLE]
for some constants , where is the wave number, is the grid size.
Proof.
Consider the difference scheme (26) as
[TABLE]
with the corresponding local truncation error given by (28).
Next, the local truncation error is the amount by which the solution value to the problem (2)-(3) does not satisfy the difference scheme (99). Hence substituting the solution value in the difference scheme (99), we get
[TABLE]
[TABLE]
or
[TABLE]
Both the constants are explicitly independent of the wave number and the solution value and given by . It is also noted that . Therefore from Lemma 4.5, we have that
[TABLE]
Thus from (102) and (103), we get the required error estimate
[TABLE]
for some constants , which completes the proof. ∎
Next we consider error estimate for scheme (19) in the following Theorem.
Theorem 4.7**.**
Let , and be a solution to the problem (2)-(3) and be the solution value to the standard sixth order scheme (19). Then for sufficiently small ,
[TABLE]
where is the wave number and is the grid size.
Proof.
Consider the scheme (19) with its corresponding local truncation error (20) and then proceeding in a similar way as for Theorem 4.6, it follows that
[TABLE]
where . Then, from Lemma 4.5,
[TABLE]
where , which completes the proof. ∎
Thus, Theorems 4.6-4.7 show that the bound of the error norm is explicitly independent of the wave number for our present scheme (26). Hence, we expect that the present scheme will provide better results for very large wave number .
Now we conclude the convergence analysis of the present difference scheme (26) in the following Theorem.
Theorem 4.8**.**
Let be a rectangular domain, with its boundary . Let and be a solution to the problem (2)-(3) and be the solution value to the present scheme (26). Then for sufficiently small , the present scheme provides sixth order accuracy. Further, as mesh size such that is sufficiently small, the present scheme converges to the solution of the corresponding boundary value problem.
5 Results and Discussion
In this section, we test the accuracy of the present difference schemes (26) and (57) over the standard sixth-order schemes [6, 11] in two and three dimensions respectively. For the purpose, several model problems governed by the Helmholtz equations in two and three dimensions are presented. All computations were performed with a 3.7 GHz Intel Xeon processor with 64 GB RAM. The program is implemented in MATLAB 2019 using BiCGstab(2) iterative algorithm with some given tolerance. Norm of the residual error is used as the stopping criteria for the iteration of the algorithm.
5.1 Two-dimensional test cases
We consider here four test case studies for the two dimensional Helmholtz equation. In order to compare the performance of the scheme (26) over the scheme [6] given by (19), we computed the - and -norms of the error and the order of convergence . We also plotted numerical errors versus different wave numbers for smaller and larger grids. If the numerical solution and exact solution values to (2)-(3) are and respectively, =grid points in - and - coordinate directions. Then the error-norms are given by
[TABLE]
and the computational order of convergence
[TABLE]
We compared the results obtained by the new scheme (26) to the 6th order scheme [6] for large wave numbers .
Problem 5.1**.**
Consider the two-dimensional Helmholtz equation [20]
[TABLE]
where , with the Dirichlet boundary
[TABLE]
Where , is odd number. The analytic solution is given by
[TABLE]
The error norms for the new scheme (26) and the sixth-order [6] for and very large are given in Tables 1-2. For this test case, the tolerance for the iteration stopping criteria is taken as . In Table 3, the error ratio of the new scheme (26) to the sixth-order scheme [6] is obtained for and . In Table 4, - error for new scheme (26) and the sixth order scheme [6] for different grid points and wave numbers is also computed.
It is noted that the source function in Problem 5.1 is independent of the wave number . From Table 1-2, it is clear that the result obtained by the new scheme is more accurate than the standard sixth order scheme. Further, Table 2 shows that the new scheme is highly accurate for large enough wave number . Table 3 shows that the given error ratio also decreases for large wave number which shows that the new scheme(26) is efficient than the scheme [6] for large wave number and small step size so that is sufficiently small.
In Figure 2, numerical error and order of accuracy of both the schemes for different wave numbers is plotted. From Figure 2(a), It is clear that the numerical errors are growing with but the numerical errors developed by the new scheme is less than that of the standard scheme. From Figure 2(b), it is also clear that the slope of each line is six in both cases for small wave number and in case of large , the slope of the line for the new scheme is more than that of the scheme [6] which implies that the scheme (26) has more accuracy than that of scheme [6] in case of large enough wave number, which also confirms the accuracy of the new scheme. The exact and numerical solutions are also plotted in Figure 3 for this test case.
In Table 4, the pollution formula is verified for the Problem 5.1, where is the order of the finite difference scheme (p=6 in this case), N is the no of mesh grids along - direction and is a constant that depends on the accuracy of the scheme. Table 4 shows that the number of grid points required for a given accuracy increases with the wave number [14, 16]. The constant is calculated by the base value and . From the Table 4, it is clear that the new scheme (26) provides high accuracy with the validation of the pollution formula.
Problem 5.2**.**
We consider the test case in Problem 5.1 with the Neumann boundary condition
[TABLE]
with the same analytic solution as in Problem 5.1, where is an odd number such that .
In order to get the sixth order accuracy, the sixth order approximation (62) for the Neumann boundary is used with the difference schemes. The maximum error norms and computational order of convergence for the Neumann problem 5.2 are computed in Tables 5-6. For this test case, the tolerance for the iteration stopping criteria is taken as . Table 5 shows that both the schemes are sixth order accurate and the new scheme is more accurate than (19). From Table 6, it is also clear that the new scheme (26) is highly accurate in case of large wave numbers with sufficiently small . We also plotted the exact and numerical solution values in Figure 4 for the Problem 5.2.
Problem 5.3**.**
Consider the model problem [20]
[TABLE]
where is a square domain with boundary . The analytic solution is given by .
The maximum error norms and computational order of convergence for the scheme (26) and (19) are computed in Tables 7-8. Table 7 shows that both the schemes are sixth order accurate and the new scheme is comparatively more accurate than (19). From Table 8, it is also clear that the new scheme (26) maintained its sixth order convergence rate even for large enough wave numbers such that is sufficiently small. From Table 7, it is also observed that the computed result obtained by both schemes are much more accurate than those obtained by fourth order difference scheme [20].
We plotted the accuracy and numerical errors for the new scheme and the standard sixth order scheme [6] in Figures 6-7. From Figure 6, it is clear that the slope of each line is six. In Figure 7, it is shown that the numerical errors are growing with the increasing for both the scheme and the numerical errors developed by the new scheme is less than that of the scheme [6], which confirms the accuracy of the new scheme. We also plotted the exact and numerical solution values in Figure 5.
Problem 5.4**.**
Consider the model problem [17]
[TABLE]
with the boundary conditions . Where are positive integers. The analytic solution is given by .
The error norms and computational order of convergence for the scheme (26) and scheme [6] are computed in Tables 9-10. Table 9 shows that both the schemes are sixth order accurate and the new scheme is more accurate than (19). From Table 10, it is also clear that the new scheme (26) is highly accurate for large wave numbers such that is sufficiently small.
The numerical errors and accuracy of the new scheme are plotted in Figures 8-9. Figure 8, clearly shows that the numerical errors are growing with the increasing for both the scheme and the numerical errors developed by the new scheme is less than that of the scheme [6]. From Figure 9, it is clear that the slope of the line for new scheme is more than that of the scheme [6] for large wave number which confirms the accuracy of the new scheme for the Problem 5.4 in case of large . The solution values for both the exact and numerical solution values are also plotted in Figure 10 for the Problem 5.4.
5.2 Three-dimensional test cases
For the three dimensional case, three test cases are presented for the Helmholtz equation in three dimensions. Hybrid BiCG method is used to solve system of equations. Results obtained from the new scheme (57) are compared with the sixth order scheme [11] given by (48) using some error norms and order of convergence. The formulae (108)-(109) can be extended in three dimensions to calculate the error norms and order of convergence for the difference schemes. Graphically, it is also shown that the numerical errors are developed with increasing the wave numbers . The exact and numerical solutions for the corresponding difference schemes are also plotted at the given grids.
Problem 5.5**.**
We extended the test problem from Fu[20] into three dimensions with the analytic solution as
[TABLE]
and the source function as
[TABLE]
where , , is odd number. The Dirichlet boundary conditions can be obtained from the analytic solution.
-error and -error norms for the new scheme (57) and the scheme (48) are given in Tables 11-12 for wave numbers . In Figure 11, numerical errors are also plotted for both the scheme with small and very large wave numbers. For this test case, iteration stops when the residual error norm is reduced by .
From tables and figures, it is noted that the numerical errors are developed for both the schemes however new scheme (57) shows more accuracy than that of the scheme (48). Further, it is also noted that the source function, in this test case, does not depend upon the wave number and therefore from Table 12, it is observed that the new scheme is highly accurate for very large wave number. The solution values for both the exact and numerical solution values are plotted in Figure 12 for the Problem 5.5.
Problem 5.6**.**
Consider the Neumann problem from [6], which is extended into three dimensions as follows
[TABLE]
with the Neumann boundary conditions
[TABLE]
where is a cubic domain. The analytic solution for the problem is given by
[TABLE]
In order to get the sixth order accuracy, it is used sixth order approximation (66) for the Neumann boundary. Maximum error norms for both the schemes (57) and (48) are given in Table 13 for different grids. For this test case, the tolerance for the iteration stopping criteria is taken . Table shows that the scheme (57) maintains its sixth order accuracy and it is also more accurate than (48). The exact and numerical solution values for the mesh grids and wave number are plotted in Figure 13 for the Problem 5.6.
Problem 5.7**.**
We consider third 3D test problem with the following analytical solution for the acoustical field [30]
[TABLE]
for which the source function is given as
[TABLE]
The Dirichlet Boundary conditions on all sides of the cubic domain can be calculated from the analytic solution.
Table 14 contains -error and maximum error norms for the schemes (57) and (48) for the Problem 5.7 with the wave number , and parameter values . Table 14 justifies the accuracy of the new scheme in comparison to the scheme (48). The tolerance for the residual error norm is taken for this test case. The exact and numerical solution values are also plotted in Figure 14 for the Problem 5.7.
6 Conclusions
A new compact sixth-order accurate finite difference scheme for the two and three-dimensional Helmholtz equation is presented. The leading truncation error term of the new scheme does not explicitly depend on the wave number. Thus the new scheme also works for problems with very large wave number . It is also shown that the new scheme is uniquely solvable for sufficiently small . Theoretically, it is proved that the bound of the error norm is explicitly independent of the wave number for the new scheme. Required symbolic derivation is performed using MAPLE with the help of mtaylor command for multi-dimensional Taylor series expansion. The resulting system of equations obtained from difference scheme are solved using BiCGstab(2) iterative method. The new scheme is tested to some model problems governed by the two and three-dimensional Helmholtz equations. Comparison of the new scheme is done with the standard sixth-order schemes [6, 11]. From the results it is shown that the new scheme is highly accurate for very high wave number. This approach can be extended to derive high-order difference schemes for the two and three dimensional problems with variable wave numbers.
Acknowledgements
The research work to the first author is supported by SRM Institute of Science and Technology Kattankulathur, Tamil Nadu India. The authors also acknowledge the SERB, New Delhi, India for the financial support towards computational facility through project file No. MTR/2017/000187.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Qing Fang, Takuya Tsuchiya, and Tetsuro Yamamoto. Finite difference, finite element and finite volume methods applied to two-point boundary value problems. Journal of Computational and Applied Mathematics , 139(1):9–19, 2002.
- 2[2] Joaquim Peiró and Spencer Sherwin. Finite difference, finite element and finite volume methods for partial differential equations. In Handbook of materials modeling , pages 2415–2446. Springer, 2005.
- 3[3] William F Ames. Numerical methods for partial differential equations . Academic press, 2014.
- 4[4] Leslaw K Bieniasz. Two new compact finite-difference schemes for the solution of boundary value problems in second-order non-linear ordinary differential equations, using non-uniform grids. Journal of Computational Methods in Sciences and Engineering , 8(1, 2):3–18, 2008.
- 5[5] M. M. Chawla. A sixth-order tridiagonal finite difference method for general non-linear two-point boundary value problems. IMA Journal of Applied Mathematics , 24(1):35–42, 1979.
- 6[6] Majid Nabavi, M. H. Kamran Siddiqui, and Javad Dargahi. A new 9-point sixth-order accurate compact finite-difference method for the helmholtz equation. Journal of Sound and Vibration , 307(3):972–982, 2007.
- 7[7] Isaac Harari and Eli Turkel. Accurate finite difference methods for time-harmonic wave propagation. Journal of Computational Physics , 119(2):252–270, 1995.
- 8[8] RP Shaw. Integral equation methods in acoustics. Boundary elements X , 4:221–244, 1988.