﻿ A Fast Fourth-Order Method for 3D Helmholtz Equation with Neumann Boundary

American Journal of Computational Mathematics
Vol.08 No.03(2018), Article ID:87259,11 pages
10.4236/ajcm.2018.83018

A Fast Fourth-Order Method for 3D Helmholtz Equation with Neumann Boundary

Na Zhu, Meiling Zhao

School of Mathematics and Physics, North China Electric Power University Department Name of Organization, Baoding, China

Received: July 28, 2018; Accepted: September 9, 2018; Published: September 12, 2018

ABSTRACT

We present fast fourth-order finite difference scheme for 3D Helmholtz equation with Neumann boundary condition. We employ the discrete Fourier transform operator and divide the problem into some independent subproblems. By means of the Gaussian elimination in the vertical direction, the problem is reduced into a small system on the top layer of the domain. The procedure for solving the numerical solutions is accelerated by the sparsity of Fourier operator under the space complexity of $O\left({M}^{3}\right)$ . Furthermore, the method makes it possible to solve the 3D Helmholtz equation with large grid number. The accuracy and efficiency of the method are validated by two test examples which have exact solutions.

Keywords:

Helmholtz Equation, Fourier Transform, Neumann Boundary Condition

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

$\begin{array}{l}\Delta \varphi +{k}^{2}\varphi =f,\text{in}\text{\hspace{0.17em}}\Omega \\ \varphi =b\left(x,y,z\right),\text{on}\text{\hspace{0.17em}}\partial \Omega \\Gamma \end{array}$ (1)

in the cubic domain $\Omega$ with Neumann boundary condition

$\frac{\partial \varphi }{\partial n}=g\left(x,y\right),\text{on}\text{\hspace{0.17em}}\Gamma ,$ (2)

where k is the wave number and $\Gamma$ is one of the planes of domain. $f\left(x,y,z\right),b\left(x,y,z\right)$ and $g\left(x,y\right)$ are known function. The Helmholtz equation is approximated by a fourth-order finite difference discretization with $h=\Delta x=\Delta y=\Delta z$ and the partition ${\left\{\left({x}_{i},{y}_{j},{z}_{l}\right)\right\}}_{i,j,l=0}^{M+1,N+1,L+1}$ .

The 19-points finite difference stencil with h yields the following linear system

$\begin{array}{l}\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\left({\delta }_{x}^{2}+{\delta }_{y}^{2}+{\delta }_{z}^{2}\right){\varphi }_{i,j,l}+\frac{{h}^{2}}{6}\left({\delta }_{x}^{2}{\delta }_{y}^{2}+{\delta }_{x}^{2}{\delta }_{z}^{2}+{\delta }_{y}^{2}{\delta }_{z}^{2}\right){\varphi }_{i,j,l}+{k}^{2}{\varphi }_{i,j,l}\\ ={f}_{i,j,l}+\frac{{h}^{2}}{12}\left({\delta }_{x}^{2}{\delta }_{y}^{2}+{\delta }_{x}^{2}{\delta }_{z}^{2}+{\delta }_{y}^{2}{\delta }_{z}^{2}\right){f}_{i,j,l}+O\left({h}^{4}\right)，\end{array}$ (3)

where ${\delta }_{x}^{2},{\delta }_{y}^{2}$ and ${\delta }_{z}^{2}$ are standard second order central difference operator and ${\varphi }_{i,j,l}$ is the fourth-order finite difference solution of Equation (1).

Moreover, we can write Equation (3) in the matrix form

$\begin{array}{l}\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\left({A}_{M}\otimes {I}_{N}\otimes {I}_{L}+{I}_{M}\otimes {A}_{N}\otimes {I}_{L}+{I}_{M}\otimes {I}_{N}\otimes {A}_{L}\right)\Phi \\ +\frac{{h}^{2}}{6}\left({A}_{M}\otimes {A}_{N}\otimes {I}_{L}+{I}_{M}\otimes {A}_{N}\otimes {A}_{L}+{A}_{M}\otimes {I}_{N}\otimes {A}_{L}\right)\Phi +{k}^{2}\Phi +{\Phi }_{B}\\ =F+\frac{{h}^{2}}{12}\left({A}_{M}\otimes {I}_{N}\otimes {I}_{L}+{I}_{M}\otimes {A}_{N}\otimes {I}_{L}+{I}_{M}\otimes {I}_{N}\otimes {A}_{L}\right)F+{F}_{B},\end{array}$ (4)

where

$\begin{array}{l}{A}_{M}=\frac{1}{{h}^{2}}\text{tridiag}\left(1,-2,1\right),{A}_{N}=\frac{1}{{h}^{2}}\text{tridiag}\left(1,-2,1\right),{A}_{L}=\frac{1}{{h}^{2}}\text{tridiag}\left(1,-2,1\right),\\ \Phi ={\left({\varphi }_{1,1,1},\cdots ,{\varphi }_{1,1,L},{\varphi }_{1,2,1},\cdots ,{\varphi }_{1,2,L},\cdots ,{\varphi }_{1,N,L},\cdots ,{\varphi }_{M,N,L}\right)}^{\text{T}},\\ F={\left({f}_{1,1,1},\cdots ,{f}_{1,1,L},{f}_{1,2,1},\cdots ,{f}_{1,2,L},\cdots ,{f}_{1,N,L},\cdots ,{f}_{M,N,L}\right)}^{\text{T}},\end{array}$

the symbol $\otimes$ represents the Kronecker product. ${I}_{M},{I}_{N},{I}_{L}$ and ${I}_{MNL}$ are identity matrices, the subscripts denote their dimension. ${A}_{M},{A}_{N}$ and ${A}_{L}$ are $M×M,N×N$ and $L×L$ tridiagonal matrices respectively. ${\Phi }_{B}$ and ${F}_{B}$ are the boundary parts of $\Phi$ and F.

3. Fast Algorithm for Three-Dimensional Helmholtz Equation

${A}_{M}$ and ${A}_{N}$ are all tridiagonal Toeplitz matrices. Fourier-sine transformation can be applied to these matrices for accelerating the algorithm. Multiplying discrete Fourier-sine transformation matrices ${S}_{M}$ and ${S}_{N}$ on the both side of ${A}_{M}$ and ${A}_{N}$ , we have

${S}_{M}{A}_{M}{S}_{M}={\text{Λ}}_{1}=\text{diag}\left({\lambda }_{1},{\lambda }_{2},\cdots ,{\lambda }_{M}\right),{S}_{N}{A}_{N}{S}_{N}={\text{Λ}}_{2}=\text{diag}\left({\mu }_{1},{\mu }_{2},\cdots ,{\mu }_{N}\right)$ ,

where

${\left({S}_{M}\right)}_{i,j}=\sqrt{\frac{2}{M+1}}\left(\mathrm{sin}\frac{ij\text{π}}{M+1}\right),{\lambda }_{i}=-\frac{4{\left(M+1\right)}^{2}}{a}{\mathrm{sin}}^{2}\frac{i\text{π}}{2\left(M+1\right)},1\le i,j\le M.$

${S}_{N}$ and ${\mu }_{t},t=1,2,\cdots ,N$ can be defined in the similar way.

Therefore, multiplying ${S}_{M}\otimes {S}_{N}\otimes {I}_{L}$ on both side of Equation (4), we have

$\begin{array}{l}\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\left({\text{Λ}}_{1}\otimes {I}_{N}\otimes {I}_{L}+{I}_{M}\otimes {\text{Λ}}_{2}\otimes {I}_{L}+{I}_{M}\otimes {I}_{N}\otimes {A}_{L}\text{}\right)\stackrel{¯}{\Phi }\\ +\frac{{h}^{2}}{6}\left({\text{Λ}}_{1}\otimes {\text{Λ}}_{2}\otimes {I}_{L}+{I}_{M}\otimes {\text{Λ}}_{2}\otimes {A}_{L}+{\text{Λ}}_{1}\otimes {I}_{N}\otimes {A}_{L}\right)\stackrel{¯}{\Phi }+{k}^{2}\stackrel{¯}{\Phi }+{\stackrel{¯}{\Phi }}_{B}\\ =\stackrel{¯}{F}+\frac{{h}^{2}}{12}\left({\text{Λ}}_{1}\otimes {I}_{N}\otimes {I}_{L}+{I}_{M}\otimes {\text{Λ}}_{2}\otimes {I}_{L}+{I}_{M}\otimes {I}_{N}\otimes {A}_{L}\right)\stackrel{¯}{F}+{\stackrel{¯}{F}}_{B},\end{array}$ (5)

where

$\begin{array}{l}\stackrel{¯}{\Phi }=\left({S}_{M}\otimes {S}_{N}\otimes {I}_{L}\right)\Phi ,\stackrel{¯}{F}=\left({S}_{M}\otimes {S}_{N}\otimes {I}_{L}\right),\\ {\stackrel{¯}{\Phi }}_{B}=\left({S}_{M}\otimes {S}_{N}\otimes {I}_{L}\right){\Phi }_{B},{\stackrel{¯}{F}}_{B}=\left({S}_{M}\otimes {S}_{N}\otimes {I}_{L}\right){F}_{B}.\end{array}$

The sparse structure of ${S}_{M}\otimes {S}_{N}\otimes {I}_{L}$ is given in Figure 1 when

Figure 1. The sparse structure of ${S}_{M}\otimes {S}_{N}\otimes {I}_{L}$ with $M=N=K=3$ .

$M=N=K=3$ , where $nz$ 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

$\begin{array}{l}\left[\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right)+\frac{{h}^{2}}{6}\left({\lambda }_{i}{\mu }_{j}{I}_{L}+{\mu }_{j}{A}_{L}+{\lambda }_{i}{A}_{L}\right)+{k}^{2}\right]{\stackrel{¯}{\Phi }}_{i,j,:}\\ ={\stackrel{¯}{F}}_{i,j,:}+\frac{{h}^{2}}{12}\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right){\stackrel{¯}{F}}_{i,j,:}+{\stackrel{¯}{F}}_{{B}_{i,j,:}}-{\stackrel{¯}{\Phi }}_{{B}_{i,j,:}},i=1,2,\cdots ,M;j=1,2,\cdots ,N.\end{array}$ (6)

In this paper, we take $\Gamma$ 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 ${\stackrel{¯}{S}}_{{B}_{top}}$ which contains the parts of ${\stackrel{¯}{\varphi }}_{i,j,L+1}$ from ${\stackrel{¯}{\Phi }}_{B}$ , there follows

$\begin{array}{l}{P}_{ij}{\stackrel{¯}{\Phi }}_{i,j,:}+\left({p}_{1}+{p}_{2}{\lambda }_{i}+{p}_{2}{\mu }_{j}\right){a}_{L2}{\stackrel{¯}{\Phi }}_{i,j,L+1}\\ ={\stackrel{¯}{F}}_{i,j,:}+\frac{{h}^{2}}{12}\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right){\stackrel{¯}{F}}_{i,j,:}+{\stackrel{¯}{F}}_{{B}_{i,j,:}}-{\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}^{\left(1\right)},\end{array}$ (7)

where

$\begin{array}{l}{P}_{ij}=\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right)+\frac{{h}^{2}}{6}\left({\lambda }_{i}{\mu }_{j}{I}_{L}+{\mu }_{j}{A}_{L}+{\lambda }_{i}{A}_{L}\right)+{k}^{2}{I}_{L},\\ {\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}^{\left(1\right)}={\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}-{\stackrel{¯}{S}}_{{B}_{top}},{p}_{1}=1+\frac{{k}^{2}{h}^{2}}{12},{p}_{2}=\frac{{h}^{2}}{6},{a}_{L2}=\frac{1}{{h}^{2}}{\left(0,0,\cdots ,1\right)}^{\text{T}}.\end{array}$

Next, we use the Gaussian elimination with a row partial pivoting to solve Equation (7).

First of all, constructing a LU-decomposition for ${P}_{ij}$ , i.e. ${P}_{ij}={L}_{ij}{U}_{ij}$ , we have

$\begin{array}{l}{L}_{ij}{U}_{ij}{\stackrel{¯}{\Phi }}_{i,j,:}+\left({p}_{2}{\lambda }_{i}+{p}_{2}{\mu }_{j}+{p}_{1}\right){a}_{L2}{\stackrel{¯}{\varphi }}_{i,j,L+1}\\ ={F}_{i,j,:}+\frac{{h}^{2}}{12}\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right){F}_{i,j,:}+{\stackrel{¯}{F}}_{{B}_{i,j,:}}-{\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}^{\left(1\right)}.\end{array}$ (8)

Since ${L}_{ij}^{-1}$ is nonsingular, multiplying ${L}_{ij}^{-1}$ on both side of Equation (8), we can obtain

$\begin{array}{l}{U}_{ij}{\stackrel{¯}{\Phi }}_{i,j,:}+\left({p}_{2}{\lambda }_{i}+{p}_{2}{\mu }_{j}+{p}_{1}\right){L}_{ij}^{-1}{a}_{L2}{\stackrel{¯}{\varphi }}_{i,j,L+1}\\ ={L}_{ij}^{-1}\left({\stackrel{¯}{F}}_{i,j,:}+\frac{{h}^{2}}{12}\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right){\stackrel{¯}{F}}_{i,j,:}+{\stackrel{¯}{F}}_{{B}_{i,j,:}}-{\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}^{\left(1\right)}\right).\end{array}$ (9)

Consequently, the last equation of Equation (9) can be derived

${\alpha }_{ij}{\stackrel{¯}{\varphi }}_{i,j,L}+{\beta }_{ij}{\stackrel{¯}{\varphi }}_{i,j,L+1}={r}_{i,j,L},i=1,2,\cdots ,M;j=1,2,\cdots ,N,$ (10)

where ${\alpha }_{ij}$ is the last element of ${U}_{ij}$ , ${\beta }_{ij}$ is the last element of $\left({p}_{2}{\lambda }_{i}+{p}_{2}{\mu }_{j}+{p}_{1}\right)\cdot {L}_{ij}{a}_{L2}$ , and ${r}_{i,j,L}$ is the last element of ${L}_{ij}^{-1}{\stackrel{¯}{F}}_{i,j,:}+\frac{{h}^{2}}{12}\left({\lambda }_{i}{I}_{L}+{\mu }_{j}{I}_{L}+{A}_{L}\right){\stackrel{¯}{F}}_{i,j,:}+{\stackrel{¯}{F}}_{{B}_{i,j,:}}-{\stackrel{¯}{\Phi }}_{{B}_{i,j,:}}^{\left(1\right)}$ . Combining $M×N$ equations analogously to Equation (10), we have

${D}_{\alpha }{\stackrel{¯}{\Phi }}_{:,:,L}+{D}_{\beta }{\stackrel{¯}{\Phi }}_{:,:,L+1}={R}_{1},$ (11)

where

$\begin{array}{l}{D}_{\alpha }=\text{diag}{\left({\alpha }_{11},{\alpha }_{12},\cdots ,{\alpha }_{1N},\cdots ,{\alpha }_{M1},{\alpha }_{M2},\cdots ,{\alpha }_{MN}\right)}^{\text{T}},\\ {D}_{\beta }=\text{diag}{\left({\beta }_{11},{\beta }_{12},\cdots ,{\beta }_{1N},\cdots ,{\beta }_{M1},{\beta }_{M2},\cdots ,{\beta }_{MN}\right)}^{\text{T}},\\ {R}_{1}={\left({r}_{1,1,L},{r}_{1,2,L},\cdots ,{r}_{1,N,L},{r}_{2,1,L},{r}_{2,2,L},\cdots ,{r}_{2,N,L},\cdots ,{r}_{M,1,L},{r}_{M,2,L},\cdots ,{r}_{M,N,L}\right)}^{\text{T}}.\end{array}$

4. Discretization of Neumann Boundary Condition

The fourth-order finite difference discretization of Equation (2) can be expressed as

$\frac{\partial \varphi }{\partial n}=\frac{{\varphi }_{i,j,L+2}-{\varphi }_{i,j,L}}{2h}-\frac{{h}^{2}}{6}{\left({\varphi }_{zzz}\right)}_{i,j,L}+O\left({h}^{4}\right).$

Using the fourth-order substitution of ${\varphi }_{zzz}$ we can derive

$\begin{array}{l}\left(1+\frac{{k}^{2}{h}^{2}}{6}+\frac{{h}^{2}}{6}{\delta }_{x}^{2}+\frac{{h}^{2}}{6}{\delta }_{y}^{2}\right){\varphi }_{i,j,L+2}-\left(1+\frac{{k}^{2}{h}^{2}}{6}+\frac{{h}^{2}}{6}{\delta }_{x}^{2}+\frac{{h}^{2}}{6}{\delta }_{y}^{2}\right){\varphi }_{i,j,L}\\ =2h{g}_{ij}+\frac{{h}^{3}}{3}{\left({f}_{z}\right)}_{i,j,L+1},i=1,2,\cdots ,M;j=1,2,\cdots ,N,\end{array}$

or the matrix form

$\begin{array}{l}\left[\left(1+\frac{{k}^{2}{h}^{2}}{6}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({A}_{M}\otimes {I}_{N}\right)+\frac{{h}^{2}}{6}\left({I}_{M}\otimes {A}_{N}\right)\right]\left({\Phi }_{:,:,L+2}-{\Phi }_{:,:,L}\right)+{\Phi }_{B}^{\left(2\right)}\\ =2hg+\frac{{h}^{3}}{3}{\left({f}_{z}\right)}_{:,:,L+1}.\end{array}$ (12)

where

$\begin{array}{l}{\Phi }_{B}^{\left(2\right)}=\frac{{h}^{2}}{6}{\left({b}_{0,1,L},{b}_{0,2,L},\cdots ,{b}_{0,N,L},0,0,\cdots ,0,\cdots ,{b}_{M+1,1,L},{b}_{M+1,2,L},\cdots ,{b}_{M+1,N,L}\right)}^{\text{T}}\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\frac{{h}^{2}}{6}{\left({b}_{1,0,L},0,\cdots ,{b}_{1,N+1,L},{b}_{2,0,L},0,\cdots ,{b}_{2,N+1,L},\cdots ,{b}_{M,0,L},0,\cdots ,{b}_{M,N+1,L}\right)}^{\text{T}}\end{array}$

and ${b}_{j,j,L}=b\left({x}_{i},{y}_{j},{z}_{L}\right)$ .

Multiplying ${S}_{M}\otimes {S}_{N}$ on both side of Equation (12), we can obtain

$\left[\left(1+\frac{{k}^{2}{h}^{2}}{6}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({\Lambda }_{1}\otimes {I}_{N}\right)+\frac{{h}^{2}}{6}\left({I}_{M}\otimes {\text{Λ}}_{2}\right)\right]\left({\stackrel{¯}{\Phi }}_{:,:,L+2}-{\stackrel{¯}{\Phi }}_{:,:,L}\right)={R}_{2}-{\stackrel{¯}{\Phi }}_{B}^{\left(2\right)},$ (13)

where

${R}_{2}=\left({S}_{M}\otimes {S}_{N}\right)\left(2hg+\frac{{h}^{3}}{3}{\left({f}_{z}\right)}_{:,:,L+1}\right),{\stackrel{¯}{\Phi }}_{B}^{\left(2\right)}=\left({S}_{M}\otimes {S}_{N}\right){\Phi }_{B}^{\left(2\right)}$ .

Moreover, replacing l with $L+1$ in Equation (3), we have

$\begin{array}{l}\left[\left(\frac{{k}^{2}{h}^{4}}{12}+\frac{2{h}^{2}}{3}\right)\left({\delta }_{x}^{2}+{\delta }_{y}^{2}\right)+\frac{{h}^{4}}{6}{\delta }_{x}^{2}{\delta }_{y}^{2}+{k}^{2}{h}^{2}-2\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)\right]{\varphi }_{i,j,L+1}\\ +\left[\left(1+\frac{{k}^{2}{h}^{2}}{12}\right)+\frac{{h}^{2}}{6}\left({\delta }_{x}^{2}+{\delta }_{y}^{2}\right)\right]\left({\varphi }_{i,j,L+2}+{\varphi }_{i,j,L}\right)\\ ={h}^{2}\left[{f}_{i,j,L+1}+\frac{{h}^{2}}{12}\left({\delta }_{x}^{2}{\delta }_{y}^{2}+{\delta }_{x}^{2}{\delta }_{z}^{2}+{\delta }_{y}^{2}{\delta }_{z}^{2}\right){f}_{i,j,L+1}\right],\end{array}$ (14)

and the matrix form

$\begin{array}{l}\left[\left(\frac{{k}^{2}{h}^{4}}{12}+\frac{2{h}^{2}}{3}\right)\left({A}_{M}\otimes {I}_{N}+{I}_{M}\otimes {A}_{N}\right)+\frac{{h}^{4}}{6}\left({A}_{M}\otimes {A}_{N}\right)-\left(2-\frac{5{k}^{2}{h}^{2}}{6}\right){I}_{MN}\right]{\Phi }_{:,:,L+1}\\ +\left[\left(1+\frac{{k}^{2}{h}^{2}}{12}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({A}_{M}\otimes {I}_{N}+{I}_{M}\otimes {A}_{N}\right)\right]\left({\Phi }_{:,:,L+2}+{\Phi }_{:,:,L}\right)+{\Phi }_{B}^{\left(3\right)}\\ ={h}^{2}\left({F}^{\left(3\right)}+{F}_{B}^{\left(3\right)}\right),\end{array}$ (15)

where

$\begin{array}{l}{F}^{\left(3\right)}=\frac{{h}^{2}}{12}\left({A}_{M}\otimes {A}_{N}+{I}_{M}\otimes {A}_{N}+{A}_{M}\otimes {I}_{N}+\frac{12}{{h}^{2}}{I}_{MN}\right){F}_{:,:,L+1}\\ \text{}-\frac{{h}^{2}}{6}\left({I}_{M}\otimes {A}_{N}+{A}_{M}\otimes {I}_{N}\right){F}_{:,:,L}\\ \text{}-\frac{{h}^{2}}{6}\left({I}_{M}\otimes {A}_{N}+{A}_{M}\otimes {I}_{N}\right){F}_{:,:,L+2}.\end{array}$

Multiplying ${S}_{M}\otimes {S}_{N}$ on both side of Equation (15), there follows

$\begin{array}{l}\left[\left(\frac{{k}^{2}{h}^{4}}{12}+\frac{2{h}^{2}}{3}\right)\left({\text{Λ}}_{1}\otimes {I}_{N}+{I}_{M}\otimes {\text{Λ}}_{2}\right)+\frac{{h}^{4}}{6}\left({\text{Λ}}_{1}\otimes {\text{Λ}}_{2}\right)-\left(2-\frac{5{k}^{2}{h}^{2}}{6}\right){I}_{MN}\right]{\stackrel{¯}{\Phi }}_{:,:,L+1}\\ +\left[\left(1+\frac{{k}^{2}{h}^{2}}{12}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({\text{Λ}}_{1}\otimes {I}_{N}+{I}_{M}\otimes {\text{Λ}}_{2}\right)\right]\left({\stackrel{¯}{\Phi }}_{:,:,L+2}+{\stackrel{¯}{\Phi }}_{:,:,L}\right)={R}_{3},\end{array}$ (16)

where ${R}_{3}={h}^{2}\left({\stackrel{¯}{F}}^{\left(3\right)}+{\stackrel{¯}{F}}_{B}^{\left(3\right)}\right)$ .

Eliminating ${\stackrel{¯}{\Phi }}_{:,:,L+2}$ from Equation (13) gives

$\begin{array}{c}C{\stackrel{¯}{\Phi }}_{:,:,L+1}+2D{\stackrel{¯}{\Phi }}_{:,:,L}={R}_{3}-D{B}^{-1}\left({R}_{2}-{\stackrel{¯}{\Phi }}_{B}{}^{\left(2\right)}\right),\end{array}$ (17)

where

$\begin{array}{l}B=\left(1+\frac{{k}^{2}{h}^{2}}{6}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({\Lambda }_{1}\otimes {I}_{N}\right)+\frac{{h}^{2}}{6}\left({I}_{M}\otimes {\text{Λ}}_{2}\right),\\ C=\left(\frac{{k}^{2}{h}^{4}}{12}+\frac{2{h}^{2}}{3}\right)\left({\text{Λ}}_{1}\otimes {I}_{N}+{I}_{M}\otimes {\text{Λ}}_{2}\right)+\frac{{h}^{4}}{6}\left({\text{Λ}}_{1}\otimes {\text{Λ}}_{2}\right)-\left(2-\frac{5{k}^{2}{h}^{2}}{6}\right){I}_{MN},\\ D=\left(1+\frac{{k}^{2}{h}^{2}}{12}\right){I}_{MN}+\frac{{h}^{2}}{6}\left({\text{Λ}}_{1}\otimes {I}_{N}+{I}_{M}\otimes {\text{Λ}}_{2}\right).\end{array}$

Combining Equation (11) and Equation (17) and derive a linear system

$\begin{array}{c}A{\stackrel{¯}{\Phi }}_{:,:,L+1}=R,\end{array}$ (18)

where

$A=C-2D{D}_{\alpha }^{-1}{D}_{\beta },R={R}_{3}-D{B}^{-1}\left({R}_{2}-{\stackrel{¯}{\Phi }}_{B}^{\left(2\right)}\right)-2D{D}_{\alpha }^{-1}{R}_{1}.$

Finally, after deriving ${\stackrel{¯}{\Phi }}_{:,:,L+1}$ , we can obtain ${\stackrel{¯}{\Phi }}_{i,j,:}$ by substituting ${\stackrel{¯}{\Phi }}_{:,:,L+1}$ in Equation (7). Multiplying ${S}_{M}\otimes {S}_{N}\otimes {I}_{L}$ , 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 $\text{Ω}=\left[0,1\right]×\left[0,1\right]×\left[0,1\right]$ .

Example 1. Consider the following problem

$\begin{array}{c}u\left(x,y,z\right)=\frac{\mathrm{sin}\left(\text{π}x\right)\mathrm{sin}\left(\text{π}y\right)}{\mathrm{sinh}\left(\sqrt{2}\text{π}\right)}\left[2\mathrm{sinh}\left(\sqrt{2}\text{π}z\right)+\mathrm{sinh}\left(\sqrt{2}\text{π}\left(1-z\right)\right)\right]\end{array}$ (19)

with

$\left\{\begin{array}{l}u\left(x,y,z\right)=\mathrm{sin}\left(\text{π}x\right)\mathrm{sin}\left(\text{π}y\right),z=0,\\ u\left(x,y,z\right)=2\mathrm{sin}\left(\text{π}x\right)\mathrm{sin}\left(\text{π}y\right),z=1,\\ u\left(x,y,z\right)=0,\text{}y\in \left\{0,1\right\},z=0,\end{array}$

$f=0$ 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 $O\left({M}^{3}\right)$ , 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 ${S}_{M}\otimes {S}_{N}\otimes {S}_{L}$ and ${S}_{M}\otimes {S}_{N}\otimes {I}_{L}$ 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 $z=1}{2}$ in Figure 2 and Figure 3 respectively.

Example 2.

$\begin{array}{l}u+{k}^{2}u=\left(-2{\text{π}}^{2}\right)\mathrm{sin}\left(\text{π}x\right)\mathrm{sin}\left(\text{π}y\right)\mathrm{sin}\left(kz\right),\text{\hspace{0.17em}}\text{in}\Omega ,\\ u=0,\text{on}\partial \Omega \\Gamma ,\end{array}$

with the exact solution

$u=\mathrm{sin}\left(\text{π}x\right)\mathrm{sin}\left(\text{π}y\right)\mathrm{sin}\left(kz\right).$ (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 $M=512$ .

Figure 3. The numerical solutions of Equation (19) on the face $z=1/2$ with $M=512$ .

Figure 4. The numerical solutions of Equation (20) with $k=3\text{π}$ (left) and $k=5\text{π}$ (right).

Figure 5. The numerical solutions of Equation (20) with $k=7\text{π}$ (left) and $k=15\text{π}$ (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).

Conflicts of Interest

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

Cite this paper

Zhu, N. and Zhao, M.L. (2018) A Fast Fourth-Order Method for 3D Helmholtz Equation with Neumann Boundary. American Journal of Computational Mathematics, 8, 222-232. https://doi.org/10.4236/ajcm.2018.83018

References

1. 1. Gupta, N., Saxena, S., Rathod, V.C. and Aggarwal, P. (2011) Unicystic Ameloblastoma of the Mandible. Journal of Oral Maxillofacial Pathology, 15, 228-231. http://dx.doi.org/10.4103/0973-029X.84511

2. 2. Menditti, D., Laino, L., Marco, G.D., Rosa, A.D., Mellone, P. and Baldi, A. (2011) Unicystic Ameloblastoma of the Mandible. In Vivo, 25, 125-128.

3. 3. Kalaskar, R., Unawane, A.S., Kalaskar, A.R. and Pandilwar, P. (2011) Conservative Management of Unicystic Ameloblastoma in a Young Child: Report of Two Cases. Contemporary Clinical Dentistry, 2, 359-363. http://dx.doi.org/10.4103/0976-237X.91804

4. 4. Ramesh, R.S., Manjunath, S., Ustad, T.H., Pais, S. and Shivakumar, K. (2010) Unicystic Ameloblastoma of the Mandible—An Unusual Case Report and Review of Literature. Head & Neck Oncology, 2, 1. http://www.headandneckoncology.org/content/2/1/1 http://dx.doi.org/10.1186/1758-3284-2-1

5. 5. Dolanmaz, D., Etoz, O.A., Pampu, A., Kalayci, A. and Gunhan, O. (2011) Marsupialization of Unicystic Ameloblastoma: A Conservative Approach for Aggressive Odontogenic Tumors. Indian Journal of Dental Research, 22, 709-712. http://dx.doi.org/10.4103/0970-9290.93461

6. 6. Yildirim, G., Ataoglu, H., Kalayci, A., Ozkan, B.T., Kucuk, K. and Esen, A. (2010) Conservative Treatment Protocol for Keratocyst Odontogenic Tumour: A Follow-Up Study of 3 Cases. Journal of Oral & Maxillofacial Research, 1, e7. http://dx.doi.org/10.5037/jomr.2010.1307

7. 7. Zhou, H., Hou, R., Ma, Q., Wu, K., Ding, Y., Qui, R. and Hu, K. (2012) Secondary Healing after Removal of Large Keratocystic Odontogenic Tumor in the Mandible: Enucleation Followed by Open Packing of Iodoform Gauze. Journal of Oral and Maxillofacial Surgery, 70, 1523-1530. http://dx.doi.org/10.1016/j.joms.2011.12.021

8. 8. Scariot, R., Silva, R.V.D., Felix Jr., W.D.S., Costa, J.D. and Rebellato, N.L.B. (2012) Conservative Treatment of Ameloblastoma in Child: A Case Report. Stomatologija, Baltic Dental and Maxillofacial Journal, 14, 33-36.

9. 9. Giuliani, M., Grossi, G.B., Lajolo, C., Bisceglia, M. and Herb, K.E. (2006) Conservative Management of a Large Odontogenic Keratocyst: Report of a Case and Review of the Literature. Journal of Oral and Maxillofacial Surgery, 64, 308-316. http://dx.doi.org/10.1016/j.joms.2005.10.013

10. 10. Maria, A., Sharma, Y. and Chabbria, A. (2012) Marsupialization as a Treatment Option of a Large Odontogenic Keratocyst: A Case Report with the Review of Literature. People’s Journal of Scientific Research, 5, 46-51.

11. 11. Miloro, M., Ghali, G.E., Larsen, P.E. and Waite, P.D. (2004) Peterson’s Principles of Oral and Maxillofacial Surgery. 2nd Edition, BC Decker Inc., Hamilton, London, Vol-1: 583.

12. 12. Wood, N.K. and Goaz, P.W. (1997) Differential Diagnosis of Oral and Maxillofacial Lesions. 5th Edition, Mosby (An Imprint of Elsevier), Maryland Heights, 337-338.

13. 13. Mendenhall, W.M., Werning, J.W., Fernandes, R., Malyapa, R.S. and Mendenhall, N.P. (2007) Ameloblastoma. American Journal of Clinical Oncology, 30, 645-648. http://dx.doi.org/10.1097/COC.0b013e3181573e59

14. 14. Yunus, M., Baig, N., Haque, A.V., Aslam, A., Atique, S., Bostan, S. and Sayed, A.M. (2009) Unicystic Ameloblastoma: A Distinct Clinicopathologic Entity. Pakistan Oral and Dental Journal, 29, 9-12.

15. 15. Kim, S.G. and Jang, H.S. (2001) Ameloblastoma: A Clinical, Radiographic and Hiistopathological Analysis of 71 Cases. Oral Surgery, Oral Medicine, Oral Pathology, Oral Radiology, Endodontology, 91, 649-653.

16. 16. Hasegawa, T., Imai, Y., Takeda, D., Yasuoka, D., Ri, S., Shigeta, T., Minamikawa, T., Shibuya, Y. and Komori, T. (2013) Retrospective Study of Ameloblastomas: The Possibility of Conservative Treatment. Kobe Journal of Medical Sciences, 59, E112-E121.

17. 17. Kahairi, A., Ahmad, R.L., Islah, L.W. and Norra, H. (2008) Management of Large Mandibular Ameloblastomas—A Case Report and Literature Reviews. Archives of Orofacial Sciences, 3, 52-55.

18. 18. Mizokami, F., Murasawa, Y., Furuta, K. and Isogai, Z. (2012) Iodoform Gauze Removes Necrotic Tissues from Pressure Ulcer Wounds by Fibrinolytic Activity. Biological and Pharmaceutical Bulletin, 35, 1048-1053. http://dx.doi.org/10.1248/bpb.b11-00016

19. 19. Freedman, M. and Stassen, L.F.A. (2013) Commonly Used Topical Oral Wound Dressing Materials in Dental and Surgical Practice—A Literature Review. Journal of the Irish Dental Association, 59, 190-195.

20. 20. Hadziabdic, N., Sulejmanagic, H., Selimovic, E. and Sulejmanagic, N. (2011) Therapeutic Approach to Large Jaw Cysts. HealthMED, 5, 1793-1799.

21. 21. Bhutia, O., Choudhury, A.R., Arora, A. and Mallick, S. (2013) Management of Unicystic Ameloblastoma of the Mandible in a 5-Year Old Child. National Journal of Maxillofacial Surgery, 4, 232-234. http://dx.doi.org/10.4103/0975-5950.127658

22. 22. Singer, I. and Turkel, E. (1998) High-Order Finite Difference Methods for the Helmholtz Equation. Computer Methods in Applied Mechanics & Engineering, 163, 343-358. https://doi.org/10.1016/S0045-7825(98)00023-1

23. 23. Harari, I. and Hughes, T.J.R. (1991) Finite Element Methods for the Helmholtz Equation in an Exterior Domain: Model Problems. Computer Methods in Applied Mechanics & Engineering, 87, 59-96. https://doi.org/10.1016/0045-7825(91)90146-W

24. 24. Jin, J.M., Liu, J., Lou, Z. and Liang, C.S.T. (2003) A Fully High-Order Finite-Element Simulation of Scattering by Deep Cavities. Antennas & Propagation IEEE Transactions on, 51, 2420-2429. https://doi.org/10.1109/TAP.2003.816354

25. 25. Jin, J.M. and Liu, J. (2000) A Special Higher Order Finite-Element Method for Scattering by Deep Cavities. IEEE Transactions on Antennas & Propagation, 48, 694-703. https://doi.org/10.1109/8.855487

26. 26. Braverman, E., Israeli, M. and Averbuch, A. (1999) A Fast Spectral Solver for a 3D Helmholtz Equation. Society for Industrial and Applied Mathematics, 20, 2237-2260. https://doi.org/10.1137/S1064827598334241

27. 27. Nabavi, M., Siddiqui, M.H.K. and Dargahi, J. (2007) A New 9-Point Sixth-Order Accurate Compact Finite-Difference Method for the Helmholtz Equation. Journal of Sound & Vibration, 307, 972-982. https://doi.org/10.1016/j.jsv.2007.06.070

28. 28. Hong, P.L., Minh, T.L. and Hoang, Q.P. (2018) On a Three Dimensional Cauchy Problem for Inhomogeneous Helmholtz Equation Associated with Perturbed Wave Number. Journal of Computational and Applied Mathematics, 335, 86-98. https://doi.org/10.1016/j.cam.2017.11.042

29. 29. Hua, Q.S., Gu, Y., Qu, W.Z., Chen, W. and Zhang, C.Z. (2017) A Meshless Generalized Finite Difference Method for Inverse Cauchy Problems Associated with Three-Dimensional Inhomogeneous Helmholtz-Type Equations. Engineering Analysis with Boundary Elements, 82, 162-171. https://doi.org/10.1016/j.enganabound.2017.06.005

30. 30. Kashirin, A.A., Smagin, S.I. and Taltykina, M.Yu. (2016) Mosaic-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

31. 31. Britt, S., Tsynkov, S. and Turkel, E. (2010) A Compact Fourth Order Scheme for the Helmholtz Equation in Polar Coordinates. Journal of Scientific Computing, 45, 26-47. https://doi.org/10.1007/s10915-010-9348-3

32. 32. Ortega, G., García, I. and Garzón, G.E.M. (2013) European Conference on Parallel Processing: A Hybrid Approach for Solving the 3D Helmholtz Equation on Heterogeneous Platforms. Springer, Berlin, 8374, 198-207.

33. 33. Sutmann, G. (2007) Compact Finite Difference Schemes of Sixth Order for the Helmholtz Equation. Journal of Computational & Applied Mathematics, 203, 15-31. https://doi.org/10.1016/j.cam.2006.03.008

34. 34. Shaw, R.P. (2000) Integral Equation Methods in Acoustics. Boundary Elements, 4, 221-244.

35. 35. Poulson, J., Engquist, B., Li, S.W. and Ying, L.X. (2013) A Parallel Sweeping Preconditioner for Heterogeneous 3D Helmholtz Equations. SIAM Journal on Scientific Computing, 35, 194-212. https://doi.org/10.1137/120871985

36. 36. Singer, I. and Turkel, E. (2006) Sixth-Order Accurate Finite Difference Schemes for the Helmholtz Equation. Journal of Computational Acoustics, 14, 339-351. https://doi.org/10.1142/S0218396X06003050

37. 37. Boisvert, R.F. (1985) A Fourth-Order-Accurate Fourier Method for the Helmholtz Equation in Three Dimensions. ACM Transactions on Mathematical Software, 13, 221-234. https://doi.org/10.1145/29380.29863

38. 38. Li, C.L. and Zou, J.W. (2013) A Sixth-Order Fast Algorithm for the Electromagnetic Scattering from Large Open Cavities. Applied Mathematics and Computation, 219, 8656-8666. https://doi.org/10.1016/j.amc.2013.02.014

39. 39. Zou, J.W. and Li, C.L. (2013) A High-Order Fast Algorithm for Three-Dimensional Helmholtz Equation. Journal of Guilin University of Electronic Technology, 33, 420-424.