A Fast Fourth-Order Method for 3D Helmholtz Equation with Neumann Boundary ()
1. Introduction
Helmholtz equation appears from general conservation laws of physics and can be interpreted as wave equations. Helmholtz equation is widely applied in the scientific and engineering design problem. Many methods have been proposed for solving the Helmholtz equations, such as finite difference method [1] , finite element method [2] [3] [4] , spectral method [5] [6] and other methods [7] [8] [9] . However, the computational cost of the finite element method increases greatly for large wave number problems. Additionally, boundary element method is limited to constant-coefficients problems. Finite difference schemes provide the simplest and least expensive avenue for achieving high-order accuracy. Some high order algorithms are proposed in [10] [11] [12] [13] . In this paper, we derive a fourth-order finite difference scheme using 19 points for solving the three-dimensional Helmholtz equation.
The discretization of the fully three-dimensional Helmholtz equation contains a large number of unknowns and requires considerable memory space. The time and space complexity increase exponentially as the grid number increases. In the meantime, to maintain a given accuracy, the mesh must be refined as the wave number increases. Some parallel algorithms are presented in [14] [15] . However, this kind of parallel algorithms cannot settle the conflict between the grid number and the performance of the computer hardware.
Fast Fourier transform is a powerful technique for solving the Helmholtz equation both in two and three dimensions [16] [17] . However, fast algorithm in [18] requires much computational cost. In light of this, we propose a fast algorithm for solving the three-dimensional Helmholtz equation. The fast operator applies inexpensive transformation to break the large discretization matrix into small and independent systems. Therefore, the equation in the whole region is divided into some small equations in the vertical direction. Meanwhile, the algorithm saves much memory space and requires less computational time due to the sparsity of the fast operator. The problem is reduced on the aperture by introducing a Gaussian elimination and the Neumann boundary condition in the vertical direction.
The paper is outlined as follows. In Section 2, a fourth-order finite difference method for the Helmholtz equation is derived. In Section 3 and Section 4, a fast algorithm is proposed by the Fourier transformation and Gaussian elimination. Two numerical experiments of the fast fourth-order algorithm are presented in Section 5. The paper is concluded in Section 6.
2. Fourth-Order Finite Difference Method
The model problem is described as follows
(1)
in the cubic domain
with Neumann boundary condition
(2)
where k is the wave number and
is one of the planes of domain.
and
are known function. The Helmholtz equation is approximated by a fourth-order finite difference discretization with
and the partition
.
The 19-points finite difference stencil with h yields the following linear system
(3)
where
and
are standard second order central difference operator and
is the fourth-order finite difference solution of Equation (1).
Moreover, we can write Equation (3) in the matrix form
(4)
where
the symbol
represents the Kronecker product.
and
are identity matrices, the subscripts denote their dimension.
and
are
and
tridiagonal matrices respectively.
and
are the boundary parts of
and F.
3. Fast Algorithm for Three-Dimensional Helmholtz Equation
and
are all tridiagonal Toeplitz matrices. Fourier-sine transformation can be applied to these matrices for accelerating the algorithm. Multiplying discrete Fourier-sine transformation matrices
and
on the both side of
and
, we have
,
where
and
can be defined in the similar way.
Therefore, multiplying
on both side of Equation (4), we have
(5)
where
The sparse structure of
is given in Figure 1 when
Figure 1. The sparse structure of
with
.
, where
means the number of the unknowns. Hence, the above equation can be transformed into a block tridiagonal matrix based on the structure of the fast operator. Equation (5) can be simplified as
(6)
In this paper, we take
as the top surface of the domain and it can be extended to the general situations. Since the solutions on the other surfaces are already known, we need to extract
which contains the parts of
from
, there follows
(7)
where
Next, we use the Gaussian elimination with a row partial pivoting to solve Equation (7).
First of all, constructing a LU-decomposition for
, i.e.
, we have
(8)
Since
is nonsingular, multiplying
on both side of Equation (8), we can obtain
(9)
Consequently, the last equation of Equation (9) can be derived
(10)
where
is the last element of
,
is the last element of
, and
is the last element of
. Combining
equations analogously to Equation (10), we have
(11)
where
4. Discretization of Neumann Boundary Condition
The fourth-order finite difference discretization of Equation (2) can be expressed as
Using the fourth-order substitution of
we can derive
or the matrix form
(12)
where
and
.
Multiplying
on both side of Equation (12), we can obtain
(13)
where
.
Moreover, replacing l with
in Equation (3), we have
(14)
and the matrix form
(15)
where
Multiplying
on both side of Equation (15), there follows
(16)
where
.
Eliminating
from Equation (13) gives
(17)
where
Combining Equation (11) and Equation (17) and derive a linear system
(18)
where
Finally, after deriving
, we can obtain
by substituting
in Equation (7). Multiplying
, we can get the numerical solution of the 3D Helmholtz equation.
5. Numerical Experiments
In this section, two numerical experiments are presented to test the validity and efficiency of the proposed method. Both experiments are implemented on MATLAB. All the equations are solved by the BiCG method. Equations in the two examples are solved in a cube
.
Example 1. Consider the following problem
(19)
with
and the corresponding Neumann boundary condition can be calculated.
Table 1 fully corroborates the theoretical design rate of the convergence for the proposed method. We can see that a good accuracy (10−7) is achieved with a small number of grid points (16 - 32 in each direction). In the case of space complexity of
, the sparsity of Fourier operator accelerates the speed for solving the three-dimensional Helmholtz equation. Moreover, the comparison of the computational time of three times Fourier transformation and twice Fourier transformation are given in Table 1. Here
and
represent two different transform operators. As we can see from Table 1, the algorithm proposed in this paper saves much computational time and makes it possible to solve the equation with large grid number. Meanwhile, we give the numerical solutions of Equation (19) in the whole domain and numerical solution on the face
in Figure 2 and Figure 3 respectively.
Example 2.
with the exact solution
(20)
We give the figures of the numerical solutions U with different wave number in Figure 4 and Figure 5. As shown in Figure 4 and Figure 5, the solutions of the Helmholtz equation are highly oscillating for large wave number.
Table 1. Convergence rate and comparisons of computational time (s) for solving Example 1 with different operators.
Figure 2. The numerical solutions of Equation (19) with
.
Figure 3. The numerical solutions of Equation (19) on the face
with
.
Figure 4. The numerical solutions of Equation (20) with
(left) and
(right).
Figure 5. The numerical solutions of Equation (20) with
(left) and
(right).
6. Conclusion
We propose a fast-high order method for solving the 3D Helmholtz equation with Neumann boundary condition. Fourier operator is used to generate block-tridiagonal structure of the discretization of the Helmholtz equation. Moreover, by using the Gaussian elimination in the vertical direction, the Helmholtz equation is reduced into a linear system in the layer of the domain. The validity and efficiency of the method are tested by two numerical experiments.
Acknowledgements
This research was supported by the Nature Science Foundation of Hebei Province (No. A2016502001) and the Fundamental Research Funds for the Central Universities (No. 2018MS129).