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

Na Zhu, Meiling Zhao^*

School of Mathematics and Physics, North China Electric Power University, Baoding, China.

DOI: 10.4236/ajcm.2019.93013                

Abstract

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.

Keywords

High-Order, Helmholtz Equation, Polar Coordinates

1. Introduction

Helmholtz equation has attracted much attention in many fields such as electromagnetic cavity scattering problems [1] , wave propagation [2] and acoustic problems [3] . Many methods have been proposed to solve the Helmholtz equation in general Cartesian coordinates, such as finite difference method [4] , finite element method [5] and other methods [6] . This equation is important for both theory and applications. Its theoretical significance has a variable coefficient under the first derivative, which renders the existing fourth-order compact finite difference methods inapplicably. From the standpoint of applications, the Helmholtz equation in polar coordinates appears in many scattering problems.

In general Cartesian coordinates, high-order methods for solving the Helmholtz equation have been well developed since the accuracy is related to the amount of the grid points with the wave number increases. Manohar et al. [7] proposed second- and sixth-order finite difference schemes for solving the 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.

2. Fourth-Order Finite Difference Scheme

2.1. Fourth-Order Approximation

We consider the following two-dimensional Helmholtz equation in polar coordinates

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+k^{2}u=f,\left(r,\theta \right)\in D,$ (1)

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

$\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)=f_r=f-k^{2}u-\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2},$ (2)

$\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}=f_{\theta }=f-k^{2}u-\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right).$ (3)

Define a uniform mesh on the region $D=\left(r_a,r_b\right)\times \left({\theta }_c,{\theta }_d\right)$ . $r_i=r_a+\left(i-1\right)h_r,i=0,1,\cdots ,M+1$ , ${\theta }_j={\theta }_c+\left(j-1\right)h_{\theta },j=0,1,\cdots ,N+1$ . $u\left(r_i,{\theta }_j\right)$ denotes the fourth-order solution of u at point $\left(r_i,{\theta }_j\right)$ .

First, we employ a fourth-order approximation on the left side of Equation (2) and obtain

$\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}{\left(r\frac{\partial u}{\partial r}\right)|}_{i,j}=\frac{1}{r_i}\frac{1}{h_r}\left(r_{i+\frac{1}{2}}\frac{u_{i+1,j}-u_{i,j}}{h_r}-r_{i-\frac{1}{2}}\frac{u_{i,j}-u_{i-1,j}}{h_r}\right)\\ -\frac{h_r^2}{12}{\left(\frac{{\partial }^2{f}_r}{\partial r^2}+\frac{1}{r}\frac{\partial f_r}{\partial r}+\frac{1}{r^2}{f}_r-\frac{2}{r^3}\frac{\partial u}{\partial r}\right)|}_{i,j}+\mathcal{O}\left(h_r^4\right),\end{array}$ (4)

where

$\frac{\partial f_r}{\partial r}=\frac{\partial f}{\partial r}-k^2\frac{\partial u}{\partial r}-\frac{\partial }{\partial r}\left(\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}\right),$

$\frac{{\partial }^2{f}_r}{\partial r^2}=\frac{{\partial }^{2}f}{\partial r^2}-k^2\frac{{\partial }^{2}u}{\partial r^2}-\frac{{\partial }^2}{\partial r^2}\left(\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}\right).$

Then approximating the derivatives of u with the second-order accuracy by central differences, we have

$\begin{array}{l}{f_r|}_{i,j}=f_{i,j}-k^2{u}_{i,j}-\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}+\mathcal{O}\left(h_r^2\right),\\ {\frac{\partial f_r}{\partial r}|}_{i,j}={\frac{\partial f}{\partial r}|}_{i,j}-k^2\frac{u_{i+1,j}-u_{i-1,j}}{2h_r}-\frac{1}{2h_r}(\frac{1}{r_{i+1}^2}\frac{u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1}}{h_{\theta }^2}\\ -\frac{1}{r_{i-1}^2}\frac{u_{i-1,j+1}-2u_{i-1,j}+u_{i-1,j-1}}{h_{\theta }^2})+\mathcal{O}\left(h_r^2\right),\\ {\frac{{\partial }^2{f}_r}{\partial r^2}|}_{i,j}={\frac{{\partial }^{2}f}{\partial r^2}|}_{i,j}-k^2\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r^2}\\ -\frac{1}{h_r^2}(\frac{1}{r_{i+1}^2}\frac{u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1}}{h_{\theta }^2}-\frac{2}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\\ +\frac{1}{r_{i-1}^2}\frac{u_{i-1,j+1}-2u_{i-1,j}+u_{i-1,j-1}}{h_{\theta }^2})+\mathcal{O}\left(h_r^2\right).\end{array}$ (5)

By substituting Equation (5) into Equation (4), we can obtain a nine-point stencil approximation for $\frac{1}{r}\frac{\partial }{\partial r}{\left(r\frac{\partial u}{\partial r}\right)|}_{i,j}$ .

We now turn to the fourth-order approximation for the left side of Equation (3). Applying a standard central difference operator to the left item of Equation (3), we have

$\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}={f_{\theta }|}_{i,j},$ (6)

and the error

$\begin{array}{l}\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\\ =\frac{1}{r}\frac{{\partial }^{2}u}{\partial {\theta }^2}+\frac{1}{r}\frac{h_{\theta }^2}{12}\frac{{\partial }^{4}u}{\partial {\theta }^4}+\mathcal{O}\left(h_{\theta }^4\right).\end{array}$ (7)

Differentiating both sides of (3) with respect to $\theta$ and combing (7), we have

${\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}|}_{i,j}=\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }}-{\frac{h_{\theta }^2}{12}\frac{{\partial }^2{f}_{\theta }}{\partial {\theta }^2}|}_{i,j}+\mathcal{O}\left(h_{\theta }^4\right).$ (8)

Similarly, all derivatives of $f_{\theta }$ with respect to $\theta$ can be approximated as follows

$\begin{array}{l}{{f^{″}}_{\theta }|}_{i,j}={\frac{{\partial }^{2}f}{\partial {\theta }^2}|}_{i,j}-k^2\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\\ -\frac{1}{r_i}\frac{1}{h_{\theta }^2{h}_r^2}\{r_{i+\frac{1}{2}}\left(u_{i+1,j+1}-u_{i,j+1}\right)-r_{i-\frac{1}{2}}\left(u_{i,j+1}-u_{i-1,j+1}\right)\\ -2\left[r_{i+\frac{1}{2}}\left(u_{i+1,j}-u_{i,j}\right)-r_{i-\frac{1}{2}}\left(u_{i,j}-u_{i-1,j-1}\right)\right]\\ +r_{i+\frac{1}{2}}\left(u_{i+1,j-1}-u_{i,j-1}\right)-r_{i-\frac{1}{2}}\left(u_{i,j-1}-u_{i-1,j-1}\right)\}+\mathcal{O}\left(h_{\theta }^2\right).\end{array}$ (9)

By substituting Equation (9) into Equation (8), the fourth-order finite difference scheme for $\frac{{\partial }^{2}u}{\partial {\theta }^2}$ is obtained.

Therefore, combing Equation (4) and Equation (8), we have the overall fourth-order finite difference form for the Helmholtz equation in polar coordinates

$\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{{\partial }^{2}u}{\partial {\theta }^2}+k^{2}u\\ =\frac{1}{r_i}\frac{1}{h_r}\left(r_{i+\frac{1}{2}}\frac{u_{i+1,j}-u_{i,j}}{h_r}-r_{i-\frac{1}{2}}\frac{u_{i,j}-u_{i-1,j}}{h_r}\right)+\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\\ -\frac{h_r^2}{12}{\left(\frac{{\partial }^2{f}_r}{\partial r^2}+\frac{1}{r}\frac{\partial f_r}{\partial r}+\frac{1}{r^2}{f}_r-\frac{2}{r^3}\frac{\partial u}{\partial r}\right)|}_{i,j}-\frac{h_{\theta }^2}{12}{\frac{{\partial }^2{f}_{\theta }}{\partial {\theta }^2}|}_{i,j}+k^2{u}_{i,j}=f_{i,j}\end{array}$ (10)

Replacing the derivatives of $f_r$ and $f_{\theta }$ by Equations (5) and (9), we can derive

$\begin{array}{l}\frac{1}{r_i}\frac{1}{h_r}\left(r_{i+\frac{1}{2}}\frac{u_{i+1,j}-u_{i,j}}{h_r}-r_{i-\frac{1}{2}}\frac{u_{i,j}-u_{i-1,j}}{h_r}\right)+\frac{1}{r_i^2}\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\\ -\frac{h_r^2}{12}\left[{\frac{{\partial }^{2}f}{\partial r^2}|}_{i,j}-k^2\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_r^2}\right]\\ +\frac{1}{12h_{\theta }^2}[\frac{1}{r_{i+1}^2}\left(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1}\right)-\frac{2}{r_i^2}\left(u_{i,j+1}-2u_{i,j}+u_{i,j-1}\right)\\ +\frac{1}{r_{i-1}^2}\left(u_{i-1,j+1}-2u_{i-1,j}+u_{i-1,j-1}\right)]\end{array}$

$\begin{array}{l}-\frac{h_r^2}{12r_i}[{\frac{\partial f}{\partial r}|}_{i,j}-k^2\frac{u_{i+1,j}-u_{i-1,j}}{2h_r}-\frac{1}{2h_r{h}_{\theta }^2}(\frac{1}{r_{i+1}^2}\left(u_{i+1,j+1}-2u_{i+1,j}+u_{i+1,j-1}\right)\\ -\frac{1}{r_{i-1}^2}\left(u_{i-1,j+1}-2u_{i-1,j}+u_{i-1,j-1}\right))]\\ -\frac{h_r^2}{12r_i^2}\left(f_{i,j}-k^2{u}_{i,j}-\frac{1}{r_i^2{h}_{\theta }^2}\left(u_{i,j+1}-2u_{i,j}+u_{i,j-1}\right)\right)\\ +\frac{h_r}{12r_i^3}\left(u_{i+1,j}-u_{i-1,j}\right)-\frac{h_{\theta }^2}{12}\left({\frac{{\partial }^{2}f}{\partial {\theta }^2}|}_{i,j}-k^2\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_{\theta }^2}\right)\end{array}$

$\begin{array}{l}+\frac{1}{12h_r^2{r}_i}[r_{i+\frac{1}{2}}\left(u_{i+1,j+1}-u_{i,j+1}\right)-r_{i-\frac{1}{2}}\left(u_{i,j+1}-u_{i-1,j+1}\right)\\ -2\left(r_{i+\frac{1}{2}}\left(u_{i+1,j}-u_{i,j}\right)-r_{i-\frac{1}{2}}\left(u_{i,j}-u_{i-1,j}\right)\right)\\ +r_{i+\frac{1}{2}}\left(u_{i+1,j-1}-u_{i,j-1}\right)-r_{i-\frac{1}{2}}\left(u_{i,j-1}-u_{i-1,j-1}\right)]+k^2{u}_{i,j}=f_{i,j}.\end{array}$ (11)

For simplicity, we rewrite the above Equation (11) in the following form

$\begin{array}{l}{l}_1{u}_{i-1,j-1}+l_2{u}_{i-1,j}+l_3{u}_{i-1,j+1}+l_4{u}_{i,j-1}+l_5{u}_{i,j}+l_6{u}_{i,j+1}\\ +l_7{u}_{i+1,j-1}+l_8{u}_{i+1,j}+l_9{u}_{j+1,j+1}=b_{i,j},\end{array}$ (12)

where

$\begin{array}{l}{l}_1=l_3=\frac{1}{12h_{\theta }^2{r}_{i-1}^2}-\frac{h_r}{24h_{\theta }^2{r}_{i-1}^2}\frac{1}{r_i}+\frac{r_{i-\frac{1}{2}}}{12h_r^2}\frac{1}{r_i},\\ l_2=\frac{5r_{i-\frac{1}{2}}}{6h_r^2{r}_i}+\frac{k^2}{12}-\frac{1}{6h_{\theta }^2{r}_{i-1}^2}+\left(\frac{h_r}{12h_{\theta }^2{r}_{i-1}^2}-\frac{h_r{k}^2}{r_i}\right)\frac{1}{12}-\frac{h_r}{12r_i^3},\\ l_4=l_6=\frac{5}{6r_i^2{h}_{\theta }^2}+\frac{h_r^2}{12r_i^4{h}_{\theta }^2}-\frac{r_{i+\frac{1}{2}}+r_{i-\frac{1}{2}}}{12h_r^2{r}_i}+\frac{k^2}{12},\\ l_5=-\frac{5}{3h_{\theta }^2{r}_i^2}-\frac{5\left(r_{i+\frac{1}{2}}+r_{i-\frac{1}{2}}\right)}{6h_r^2{r}_i}+\frac{2k^2}{3}-\left(\frac{h_r^2}{6r_i^2{h}_{\theta }^2}-\frac{h_r^2{k}^2}{r_i^2}\right)\frac{1}{r_i^2},\end{array}$

$\begin{array}{l}{l}_7=l_9=\frac{1}{12h_{\theta }^2{r}_{i+1}^2}+\frac{h_r}{24h_{\theta }^2{r}_{i+1}^2}\frac{1}{r_i}+\frac{r_{i+\frac{1}{2}}}{12h_r^2}\frac{1}{r_i},\\ l_8=\frac{5r_{i+\frac{1}{2}}}{6h_r^2{r}_i}+\frac{k^2}{12}-\frac{1}{6h_{\theta }^2{r}_{i+1}^2}-\left(\frac{h_r}{12h_{\theta }^2{r}_{i+1}^2}-\frac{h_r{k}^2}{12}\right)\frac{1}{r_i}+\frac{h_r}{12r_i^3},\\ b_{i,l}=f_{i,l}+\frac{h_r^2}{12r_i^2}{f}_{i,l}+\frac{h_r^2}{12r_i}\frac{\partial f}{\partial r}+\frac{h_r^2}{12}\frac{{\partial }^{2}f}{\partial r^2}+\frac{h_{\theta }^2}{12}\frac{{\partial }^{2}f}{\partial {\theta }^2}.\end{array}$

2.2. The Sparse Matrix Form of the Scheme

The relationship between the nine points can be illustrated in three parts in matrix form $A_X,A_Y,A_D$ , which represent the horizontal, vertical, diagonal relationships respectively,

$\left(A_X+A_Y+A_D\right)U=F,$ (13)

where

$\begin{array}{l}U={\left(u_{1,1},u_{1,2},\cdots ,u_{1,N},u_{2,1},u_{2,2},\cdots ,u_{2,N},\cdots ,u_{M,1},u_{M,2},\cdots ,u_{M,N}\right)}^{\text{T}},\\ A_X=A_x\otimes I_N,A_Y=\text{diag}\left\{A_Y^{\left(1\right)},A_Y^{\left(2\right)},\cdots ,A_Y^{\left(M\right)}\right\},A_Y^{\left(i\right)}=\text{tridiag}\left\{l_6^{\left(i\right)},\frac{l_5^{\left(i\right)}}{4},l_4^{\left(i\right)}\right\},\\ A_D=A_{D_1}\otimes I_{D_1}+A_{D_2}\otimes I_{D_2}+A_{D_3}\otimes I_N,\\ A_x=\left(\begin{array}{cccc}\frac{l_5^{\left(1\right)}}{4}& l_8^{\left(1\right)}& & \\ l_2^{\left(2\right)}& \frac{l_5^{\left(2\right)}}{4}& l_8^{\left(2\right)}& \\ & & \ddots & \\ & & l_2^{\left(M\right)}& \frac{l_5^{\left(M\right)}}{4}\end{array}\right),A_{D_1}=\left(\begin{array}{ccccc}0& l_7^{\left(1\right)}& & & \\ l_1^{\left(2\right)}& 0& l_7^{\left(2\right)}& & \\ & & \ddots & \ddots & \\ & & l_1^{\left(M-1\right)}& 0& l_7^{\left(M-1\right)}\\ & & & l_1^{\left(M\right)}& 0\end{array}\right),\end{array}$

$\begin{array}{l}{I}_{r_1}=\left(\begin{array}{ccccc}0& 1& & & \\ & 0& 1& & \\ & & \ddots & \ddots & \\ & & & 0& 1\\ & & & & 0\end{array}\right),A_{D_2}=\left(\begin{array}{ccccc}0& l_9^{\left(1\right)}& & & \\ l_3^{\left(2\right)}& 0& l_9^{\left(2\right)}& & \\ & & \ddots & \ddots & \\ & & l_3^{\left(M-1\right)}& 0& l_9^{\left(M-1\right)}\\ & & & l_3^{\left(M\right)}& 0\end{array}\right),\\ I_{r_2}=\left(\begin{array}{ccccc}0& & & & \\ 1& 0& & & \\ & \ddots & \ddots & & \\ & & 1& 0& \\ & & & 1& 0\end{array}\right),A_{D_3}=\left(\begin{array}{cccc}\begin{array}{l}\frac{l_5^{\left(1\right)}}{2}\\ \\ \\ \end{array}& \begin{array}{l}\\ \frac{l_5^{\left(2\right)}}{2}\\ \\ \end{array}& \begin{array}{l}\\ \\ \ddots \\ \end{array}& \begin{array}{l}\\ \\ \\ \frac{l_5^{\left(M\right)}}{2}\end{array}\end{array}\right)\end{array}$

and $\otimes$ denotes the Kronecker product, $I_N$ is the $N\times N$ identity matrix. F is the right side item containing $b_{i,j}$ and the boundary conditions which will be discussed in the following section.

3. Boundary Condition

Implementation of Dirichlet boundary conditions at $r=r_a,r_b$ , $\theta ={\theta }_c,{\theta }_d$ is straightforward. And the right side item of Equation (13) can be written as

$F=\tilde{f}+U_{Bl}+U_{Br}+U_{Bb}+U_{Bt},$ (14)

where

$\begin{array}{l}\tilde{f}={\left(b_{11},b_{12},\cdots ,b_{1N},b_{21},b_{22},\cdots ,b_{2N},\cdots ,b_{M1},b_{M2},\cdots ,b_{MN}\right)}^{\text{T}},\\ U_{Bl}={\left(B_1{u}_{0,:},0,0,\cdots ,0,0,0,\cdots ,0,0,0,\cdots ,0\right)}^{\text{T}},u_{0,:}={\left(u_{0,0},u_{0,1},\cdots ,u_{0,N+2}\right)}^{\text{T}},\\ U_{Br}={\left(B_2{u}_{M+1,:},0,0,\cdots ,0,0,0\cdots ,0,0,0,\cdots ,0\right)}^{\text{T}},\\ u_{M,:}={\left(u_{M,0},u_{M,1},\cdots ,u_{M,N+2}\right)}^{\text{T}},U_{Bb}=B_3{u}_{:,0}\otimes a_N,U_{Bb}=B_4{u}_{:,N+1}\otimes b_N,\\ u_{:,0}={\left(u_{1,0},u_{2,0},\cdots ,u_{M,0}\right)}^{\text{T}},u_{:,N+1}={\left(u_{1,N+1},u_{2,N+1},\cdots ,u_{M,N+1}\right)}^{\text{T}},\end{array}$

$\begin{array}{l}{B}_1=\left(\begin{array}{cccccc}{l}_3^{\left(1\right)}& l_2^{\left(1\right)}& l_1^{\left(1\right)}& & & \\ & l_3^{\left(1\right)}& l_2^{\left(1\right)}& l_1^{\left(1\right)}& & \\ & & \ddots & \ddots & \ddots & \\ & & & l_3^{\left(1\right)}& l_2^{\left(1\right)}& l_1^{\left(1\right)}\end{array}\right),\\ B_2=\left(\begin{array}{cccccc}{l}_9^{\left(M\right)}& l_8^{\left(M\right)}& l_7^{\left(M\right)}& & & \\ & l_9^{\left(M\right)}& l_8^{\left(M\right)}& l_7^{\left(M\right)}& & \\ & & \ddots & \ddots & \ddots & \\ & & & l_9^{\left(M\right)}& l_8^{\left(M\right)}& l_7^{\left(M\right)}\end{array}\right),\\ B_3=B_4=\left(\begin{array}{ccccc}{l}_4^{\left(1\right)}& l_7^{\left(1\right)}& & & \\ l_1^{\left(2\right)}& l_4^{\left(2\right)}& l_7^{\left(2\right)}& & \\ & \ddots & \ddots & \ddots & \\ & & l_1^{\left(M-1\right)}& l_4^{\left(M-1\right)}& l_7^{\left(M-1\right)}\\ & & & l_4^{\left(M\right)}& l_7^{\left(M\right)}\end{array}\right),a_N=\left(\begin{array}{c}1\\ 0\\ ⋮\\ 0\end{array}\right),b_N=\left(\begin{array}{c}0\\ 0\\ ⋮\\ 1\end{array}\right),\end{array}$

and $a_N,b_N$ are vectors in $N\times 1$ dimensions.

Implementation of the Dirichlet boundary conditions at $r=r_a,r_b$ , $\theta ={\theta }_c,{\theta }_d$ is more complicated. We consider the following Neumann boundary condition and it can be extended to the general cases

${\frac{\partial u}{\partial \theta }|}_{\theta ={\theta }_{top}}=\alpha u+g\left(r,\theta \right),$ (15)

where $\alpha$ is a constant, $g\left(r,\theta \right)$ is a given function.

The ghost points are utilized to derive the difference scheme on the boundary can be written as

${\frac{\partial u}{\partial \theta }|}_{i,N+1}=\frac{u_{i,N+2}-u_{i,N}}{2h_{\theta }}+\mathcal{O}\left(h^2\right).$ (16)

Therefore, the Neumann boundary condition can be approximated as

$\alpha u_{i,N+1}+g\left(r_i,{\theta }_{N+1}\right)=\frac{u_{i,N+2}-u_{i,N}}{2h_{\theta }},i=1,2,\cdots ,M.$ (17)

We can obtain

$u_{:,N+2}-u_{:,N}=2h_{\theta }{g}_{:,N+1}+2h_{\theta }\alpha u_{:,N+1},$ (18)

where

$\begin{array}{l}{u}_{:,N}={\left(u_{1,N},u_{2,N},\cdots ,u_{M,N}\right)}^{\text{T}},\\ u_{:,N+1}={\left(u_{1,N+1},u_{2,N+1},\cdots ,u_{M,N+1}\right)}^{\text{T}},\\ u_{:,N+2}={\left(u_{1,N+2},u_{2,N+2},\cdots ,u_{M,N+2}\right)}^{\text{T}},\\ g_{:,N+1}={\left(g_{1,N+1},g_{2,N+1},\cdots ,g_{M,N+1}\right)}^{\text{T}}.\end{array}$ (19)

Substituting j for $N+1$ in Equation (12) gives

$A{u}_{:,N+2}+B{u}_{:N+1}+A{u}_{:,N}=b_{:,N+1},$ (20)

where

$\begin{array}{l}A=\left(\begin{array}{ccccc}{l}_4^{\left(1\right)}& l_7^{\left(1\right)}& & & \\ l_1^{\left(2\right)}& l_4^{\left(2\right)}& l_7^{\left(2\right)}& & \\ & \ddots & \ddots & \ddots & \\ & & l_1^{\left(M-1\right)}& l_4^{\left(M-1\right)}& l_7^{\left(M-1\right)}\\ & & & l_4^{\left(M\right)}& l_7^{\left(M\right)}\end{array}\right),\\ B=\left(\begin{array}{ccccc}{l}_5^{\left(1\right)}& l_8^{\left(1\right)}& & & \\ l_2^{\left(2\right)}& l_5^{\left(2\right)}& l_8^{\left(2\right)}& & \\ & \ddots & \ddots & \ddots & \\ & & l_2^{\left(M-1\right)}& l_5^{\left(M-1\right)}& l_8^{\left(M-1\right)}\\ & & & l_5^{\left(M\right)}& l_8^{\left(M\right)}\end{array}\right),\\ b_{:,N+1}={\left(b_{1,N+1},b_{2,N+1},\cdots ,b_{M,N+1}\right)}^{\text{T}}.\end{array}$

Combining Equations (18) and (20), we can eliminate $u_{:,N+2}$ and obtain

$2A{u}_{:,N}+\left(2h_{\theta }\alpha A+B\right)u_{:,N+1}=b_{:,N+1}-2h_{\theta }A{g}_{:,N+1}.$ (21)

Therefore, collaborating Equations (13) and (21), the global system can be written as follows

$\left(\begin{array}{cc}{A}_{11}& A_{12}\\ A_{21}& A_{22}\end{array}\right)\left(\begin{array}{c}U\\ u_{:,N+1}\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\end{array}\right),$ (22)

where

$\begin{array}{l}{A}_{11}=A_X+A_Y+A_D,A_{12}=B_4,A_{21}=2A\otimes a_N,A_{22}=2h_{\theta }\alpha A+B,\\ b_1=\tilde{f}+U_{Bl}+U_{Br}+U_{Bb},b_2=b_{:,N+1}-2h_{\theta }A{g}_{:,N+1}.\end{array}$

We can observe the sparsity of the global system of the Helmholtz equation with the Neumann boundary condition in Figure 1.

4. Numerical Experiments

4.1. Example 1

We first consider the problem with Dirichlet boundary condition in polar coordinates as follows

$\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+k^{2}u=0,\left(r,\theta \right)\in \left[1,2\right]\times \left[0,2\text{π}\right]\\ {u|}_{r=1}=\sin \left(\theta \right),{u|}_{r=2}=2\sin \left(\theta \right),{u|}_{\theta =0}={u|}_{\theta =2\text{π}}=0.\end{array}$ (23)

and the exact solution is $u\left(r,\theta \right)=r\cos \left(\theta \right)$ .

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 Table 2. The order is calculated by the following equation

$\text{order}=\log \left(\frac{\text{error}\left(M_1\right)}{\text{error}\left(M_2\right)}\right)/\log \left(\frac{M_2}{M_1}\right).$

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.

Figure 2. The numerical solution (left) and the exact solution with $k=1$ and $M=N=64$ .

Figure 3. The error inside the region (left) and the numerical solution in Cartesian coordinates (right) with $k=1$ and $M=N=64$ .

Table 1. The errors and order of the proposed method with $k=1$ and $k=5$ .

Table 2. The errors and order of the proposed method with $k=10$ and $k=20$ .

Table 3. Computational time (s) for solving the Helmholtz equation with $k=10$ .

gives the comparison of the computational time(s) for solving the Helmholtz equation in polar coordinates, where Method I and Method II denote methods with and without the utilization of the sparsity of the coefficient matrix of the linear system. It can be observed that with the help of the sparsity of the coefficient matrix of the linear system, the computational cost is remarkably reduced.

4.2. Example 2

This example is a Helmholtz equation in polar coordinates with Neumann boundary condition

$\begin{array}{l}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+k^{2}u=0,\left(r,\theta \right)\in \left[1,2\right]\times \left[0,2\text{π}\right]\\ \frac{\partial u}{\partial \theta }=r\cos \left(\theta \right),\theta =2\text{π},{u|}_{r=1}=\sin \left(\theta \right),{u|}_{r=2}=2\sin \left(\theta \right),{u|}_{\theta =0}=0,\end{array}$ (24)

and the exact solution is $u=r\sin \left(\theta \right)$ . As we can see from Figures 4-6 that the

Figure 4. The numerical solution (left) and the exact solution with $k=5$ and $M=N=128$ .

Figure 5. The error inside the region (left) and the error on the top boundary (right) with $k=5$ and $M=N=128$ .

Figure 6. The numerical solution in Cartesian coordinates.

numerical solution is well agreed with the exact solution with $k=5$ and $M=N=128$ .

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

Acknowledgements

This research was supported by the Fundamental Research Funds for the Central Universities (No. 2018MS129).

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zhao, M.L. (2013) A Fast High Order Iterative Solver for the Electromagnetic Scattering by Open Cavities Filled with Inhomogeneous Media. Advances in Applied Mathematics and Mechanics, 5, 235-257. https://doi.org/10.4208/aamm.12-m12119
[2]	Zhao, C. and Liu, T. (2003) Nonlinear Artificial Boundaries for Transient Scalar Wave Propagation in a Two-Dimensional Infinite Homogeneous Layer. International Journal for Numerical Methods in Engineering, 58, 1435-1456. https://doi.org/10.1002/nme.703
[3]	Colton, D. and Kress, R. (1998) Inverse Acoustic and Electromagnetic Scattering Theory. Springer, Berlin. https://doi.org/10.1007/978-3-662-03537-5
[4]	Singer, I. and Turkel, E. (1998) High-Order Finite Difference Methods for the Helmholtz Equation. Computer Methods in Applied Mechanics and Engineering, 163, 343-358. https://doi.org/10.1016/S0045-7825(98)00023-1
[5]	Harari, I. and Hughes, T. (1991) Finite Element Methods for the Helmholtz Equation in an Exterior Domain: Model Problems. Computer Methods in Applied Mechanics and Engineering, 87, 59-96. https://doi.org/10.1016/0045-7825(91)90146-W
[6]	Kashirin, A., Smagin, S. and Taltykina, M. (2016) Osaic-Skeleton Method as Applied to the Numerical Solution of Three-Dimensional Dirichlet Problems for the Helmholtz Equation in Integral Form. Computational Mathematics and Mathematical Physics, 56, 612-625. https://doi.org/10.1134/S0965542516040096
[7]	Manohar, R. and Tephenson, J. (1983) Single Cell High Order Difference Methods for Helmholtz Equation. Journal of Computational Physics, 51, 444-453. https://doi.org/10.1016/0021-9991(83)90163-8
[8]	Nabavi, M., Siddiqui, M. and Dargahi, J. (2007) A New 9-Point Sixth-Order Accurate Compact Finite Difference Method for the Helmholtz Equation. Journal of Sound and Vibration, 37, 972-982. https://doi.org/10.1016/j.jsv.2007.06.070
[9]	Sutmann, G. (2007) Compact Finite Difference Schemes of Sixth-Order for the Helmholtz Equation. Journal of Computational and Applied Mathematics, 203, 15-31. https://doi.org/10.1016/j.cam.2006.03.008
[10]	Zhang, R.S., Lu, G.Z. and Abomakhleb, G. (2017) Finite-Difference Solution of the Parabolic Equation under Horizontal Polar Coordinates. IEEE Antennas and Wireless Propagation Letters, 16, 2931-2934. https://doi.org/10.1109/LAWP.2017.2753220
[11]	Britt, S., Tsynkov, S. and Turkel, E. (2001) Numerical Simulation of Time-Harmonic Waves in Inhomogeneous Media Using a Compact High Order Schemes. Communications in Computational Physics, 9, 520-541. https://doi.org/10.4208/cicp.091209.080410s
[12]	Su, X.L., Feng, X.F. and Li, Z.L. (2016) Fourth-Order Compact Schemes for Helmholtz Equation with Piecewise Wave Numbers in the Polar Coordinates. Journal of Computational Mathematics, 34, 499-510. https://doi.org/10.4208/jcm.1604-m2015-0290
[13]	Jha, N., Monhanty, R.K. and Kumar, N. (2017) Compact-FDM for Midly Nonlinear Two-Space Dimensional Elliptic BVPs in Polar Coordinate System and Its Convergence Theory. Journal of Computational and Applied Mathematics, 3, 255-270. https://doi.org/10.1007/s40819-015-0104-0