A Fast High Order Algorithm for 3D Helmholtz Equation with Dirichlet Boundary

Helmholtz equation is widely applied in the scientific and engineering problem. For the solution of the threedimensional Helmholtz equation, the computational efficiency of the algorithm is especially important. In this paper, in order to solve the contradiction between accuracy and efficiency, a fast high order finite difference method is proposed for solving the three-dimensional Helmholtz equation. First, a traditional fourth order method is constructed. Then, fast Fourier transformation are introduced to generate a block-tridiagonal structure which can easily divide the original problem into small and independent subsystems. For large 3D problems, the computation of traditional discrete Fourier transformation is less efficient, and the memory requirements increase rapidly. To fix this problem, the transformed coefficient matrix is constructed as a sparse structure. In light of the sparsity, the algorithm presented in this paper requires less memory space and computational cost. This sparse structure also leads to independent solving procedure of any plane in the domain. Therefore, parallel method can be used to solve the Helmholtz equation with large grid number. In the end, three numerical experiments are presented to verify the effectiveness of the fast fourth-order algorithm, and the acceleration effect to use the parallel method has been demonstrated.


Introduction
In this paper, we consider the following Helmholtz equation with Dirichlet boundary condition , , , , , on Ω, # where is the wave number. Helmholtz equation is widely applied to wave propagation and acoustic scattering problem, viscous and heat transfer problems and electromagnetic scattering problems.
In recent years, there has been tremendous interest in developing algorithms for solving the Helmholtz equation. The precise numerical solution of Helmholtz equation is of great importance. Most algorithms were initially developed for the two-dimensional case, such as finite element method, finite difference method and other iterative methods [1][2][3][4]. However, the accuracy of the finite element method reduces rapidly with the increasement of the wave number. In addition, finite difference method is not competitive for solving the Helmholtz equation with large wave number if no fine mesh provided. Therefore, some high order finite difference method were proposed for solving the twodimensional Helmholtz equation in [5][6][7][8][9][10]. This method is prevalent since they can offer a high accuracy solutions with less computational cost. There have been applied for solving the three-dimensional Helmholtz equations. Braverman et al. proposed a fast spectral Helmholtz solver which incorporates the application of the FFT with a preliminary subtraction technique in [11]. Lu presented a fourth-order compact difference scheme based on the Laplace operator in [12]. A noniterative solver based on the domain decomposition was developed in [13] for solving the three-dimensional Helmholtz equation. Sutmann combined different high order finite difference schemes and proposed a more efficient algorithm in [14,15]. Furthermore, Gordon et al. extend the sixth-order finite difference method to high and variable wave number with CARP-CG algorithm in [16,17]. Due to the heavy computational cost the efficiency of the algorithm for solving three-dimensional Helmholtz equation is especially important. The Fourier method has many advantages for solving the Helmholtz equation. Together with the symmetry of the high-order finite difference operator, the Dirichlet problems can be solved easily with sine transformation.
In this paper, our aim is to discretize the Helmholtz equation with fourth-order finite difference method and build a fast algorithm using the discrete Fourier-sine transformation. On the other hand, the traditional Fourier transformation in three dimensions needs considerable memory requirement especially for large wave number equations. In views of this problem, twice Fourier transformation are applied for the discrete Helmholtz equation.
The paper is organized as follows. In Section 2, a fourthorder finite difference method was proposed to discrete the three-dimensional Helmholtz equation. In Sections 3 and 4, a fast high algorithm is developed for solving the threedimensional Helmholtz equation. Section 5 describes the parallel implementation for the three-dimensional Helmholtz equation. Moreover, three numerical examples prove the validity and feasibility of the proposed method. The paper is concluded in Section 6.

The Fourth-Order Finite Difference Method
We consider a fourth-order finite difference method of (1), where Ω is the model domain and Ω is the boundary of the domain.
, , is an given continuous function. In this paper, a cubic domain Ω 0, 0, 0, is considered as the solution field. Firstly, the uniform partition is define as , , ! , ," #$%,&$%,'$% in Ω . Without loss of generality, we consider the case of ( ∆ ∆ ∆ since it can be extended to the general situation. For any point , , in Ω the standard second order central difference operator can be written as where    R H XYS can be written as follows

The Fast Algorithm for 3D Helmholtz Equation
For the tridiagonal Toeplitz matrix D # and D & , we have R # and R & are discrete Fourier-sine transformation matrices. R & and k S , p = 1,2, … , 5 can be defined similarly as (8). Multiplying R # ⊗ R & ⊗ F ' on both side of (4), the following formula can be obtained

Preparation for the Numerical Experiments
In this section, the algorithm is applied for three numerical experiments. And the performance of the algorithm is described using MATLAB on an eight processor computer.
One error measurement is used to weigh the difference between , ,9 * and , ,9 , where , ,9 * denotes the real solution of the three-dimensional Helmholtz equation.
Define the error measurement as follows

Parallel Implementation
After deriving G q by (10), the solution G of the threedimensional Helmholtz equation is calculated by G R # E R & E F ' G q . Due to the special form of the coefficient matrix, the computational procedure is independent for each , 6 0,1, … , 7 . Therefore, the fast algorithm proposed in Section 3 is well adapted for the parallel computation. The implementation under … processors is shown in Figure 3.
• , , sin OE sin OE , 0 , , , 2 sin OE sin OE , 1, 0, , Ž •0,1•, The numerical results for solving G q and G are displayed in Table 1. Here CPU " q and CPU " denotes the computational time (s) for solving G q and G respectively. After deriving G q , we need to multiply R # E R & E F ' on the left side of G q to compute the numerical solution of the Helmholtz equation. As we can see from Table 1, CPU " q time is notably less than CPU " time especially for large gird points. Moreover, the comparison of the computational time of three times Fourier transformation and twice Fourier transformation are given in Table 1. Here R # E R & E R ' and R # E R & E F ' represent two different transform operators.    The behavior of the parallel implementation was examined in Example 2. The results are displayed in Table 2. The results again verify the efficiency of the parallel processor. Moreover, we test the convergence order of the proposed fast high algorithm. The last column of Table 2 demonstrates the fourth-order convergence of the proposed algorithm. The line in Figure 6 is fitted by log error and log 3 and its slope is 4.0545.
The figures of the numerical solutions G with different wave number in Figure 7 and

Conclusion
In this paper, we proposed an fourth-order fast algorithm for solving the three-dimensional Helmholtz equation. By multiplying the Fourier transform operator, the large system becomes some small independent systems. Moreover, in view of the special format of the coefficient matrix, we utilize the parallel processors to solve the equation. This implementation is especially efficient for solving the numerical solutions on one of the surfaces. Three numerical experiments have demonstrated the validity of the fourth-order fast algorithm.