High-Order Finite Difference Method for Helmholtz Equation in Polar Coordinates

We present a fourth-order finite difference scheme for the Helmholtz equation in polar coordinates. We employ the finite difference format in the interior of the region and derive a nine-point fourth-order scheme. Specially, ghost points outside the region are applied to obtain the approximation for the Neumann boundary condition. We obtain the matrix form of the linear system and the sparsity of the coefficient matrix is favorable for the computation of the Helmholtz equation. The feasibility and accuracy of the method are validated by two test examples which have exact solutions.

Helmholtz equation which requires the derivatives and much computational cost. Nabavi et al. [8] further developed a compact nine-point sixth-order finite difference scheme and obtained a sixth-order approximation for the Neumann boundary condition. Sutmann [9] developed a new compact sixth-order finite difference scheme of sixth-order for three-dimensional Helmholtz equations.
In polar coordinates, the symmetric problem can be described more concisely. A finite difference method was proposed to solve the parabolic equation in polar coordinates in [10]. Britt et al. [11] constructed a high-order compact difference scheme for the Helmholtz equation. Su et al. [12] proposed a fourth-order method for solving Helmholtz equation with discontinuous wave number and Dirichlet boundary condition. A novel compact scheme based on finite difference discretizations and geometric grid has been developed to solve two-dimensional mildly non-linear elliptic equations in polar co-ordinate constituting singular terms [13].
The existed works didn't discuss the Neumann boundary condition and the sparsity of the coefficient matrix. In this paper, we develop a fourth-order scheme for the Helmholtz equation in polar coordinates with Dirichlet and Neumann boundary conditions. The sparse matrix form of the linear system is derived. With the help of the sparsity of the coefficient matrix of the linear system, the computational cost is remarkably reduced. Moreover, we give the approximation for the Neumann boundary condition. The feasibility and the order of method are verified by the Helmholtz equation with Dirichlet and Neumann boundary condition.
The paper is outlined as follows. In Section 2, a fourth-order finite difference scheme and the sparse matrix form for the Helmholtz equation in polar coordinates are derived. In Section 3, the approximation for the boundary condition is obtained. Two numerical experiments of the high-order algorithm are presented in Section 4. The paper is concluded in Section 5.

Fourth-Order Approximation
We consider the following two-dimensional Helmholtz equation in polar coordinates ( ) where k is the wave number, D is a given region and f is a known function. For convenience, we rewrite the above equation as 1 , 1 .
Then approximating the derivatives of u with the second-order accuracy by By substituting Equation (5) into Equation (4), we can obtain a nine-point stencil approximation for We now turn to the fourth-order approximation for the left side of Equation Differentiating both sides of (3) with respect to θ and combing (7), we have ( ) Similarly, all derivatives of f θ with respect to θ can be approximated as By substituting Equation (9) into Equation (8), the fourth-order finite difference scheme for Therefore, combing Equation (4) and Equation (8), we have the overall fourth-order finite difference form for the Helmholtz equation in polar coordinates Replacing the derivatives of r f and f θ by Equations (5) and (9), we can For simplicity, we rewrite the above Equation (11) in the following form

The Sparse Matrix Form of the Scheme
The relationship between the nine points can be illustrated in three parts in ma- and ⊗ denotes the Kronecker product, N I is the N N × identity matrix. F is the right side item containing , i j b and the boundary conditions which will be discussed in the following section.

Implementation of Dirichlet boundary conditions at
where α is a constant, ( ) , g r θ is a given function.
The ghost points are utilized to derive the difference scheme on the boundary can be written as ( ) Therefore, the Neumann boundary condition can be approximated as Therefore, collaborating Equations (13) and (21), the global system can be written as follows 11 12 1

Example 1
We first consider the problem with Dirichlet boundary condition in polar coordinates as follows and the exact solution is ( ) ( ) , cos u r r θ θ = . As we can see from Figure 2 and Figure 3, the numerical solution derived by the proposed method is highly consistent with the exact solution. Moreover, we depict the numerical solution in Cartesian coordinates in the right side of Figure  3. Furthermore, in order to test the computational order of the proposed method, we give the error between the numerical solution and the exact solution with different grid points and k in Table 1 and It can be clearly seen from Table 1 and Table 2 that the finite difference scheme can reach the fourth-order when the wavenumber k is relatively small.
As the mesh is refined and the number of grid points increases, the error becomes smaller and the accuracy tends to be fourth-order and gradually stabilized. When the wavenumber k increases, the error becomes oscillatory. Table 3 Figure 1. The sparsity of the global system with Neumann boundary condition.

Conclusion
In this paper, we propose a high-order fast algorithm for solving the two-dimensional Helmholtz equation with Dirichlet and Neumann boundary conditions in polar coordinates. We develop a fourth-order accurate compact finite difference approximation to the Helmholtz equation. The sparse matrix form for the Helmholtz equation in polar coordinates is obtained which improves the efficiency for the computation process. Two numerical experiments have demonstrated the validity of the fourth-order algorithm.