# Multigrid Method for Solution of 3D Helmholtz Equation Based on HOC Schemes

Authors: Fazal Ghaffar; Noor Badshah; Saeed Islam
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A higher order compact difference (HOC) scheme with uniform mesh sizes in different coordinate directions is employed to discretize a two- and three-dimensional Helmholtz equation. In case of two dimension, the stencil is of 9 points
while in three-dimensional case, the scheme has 27 points and has fourth- to fifth-order accuracy. Multigrid method using
Gauss-Seidel relaxation is designed to solve the resulting sparse linear systems. Numerical experiments were conducted
to test the accuracy of the sixth-order compact difference scheme with Multigrid method and to compare it with the
standard second-order finite-difference scheme and fourth-order compact difference scheme. Performance
of the scheme is tested through numerical examples. Accuracy and efficiency of the new scheme are established
by using the errors normsl
2.


---

## Body

## 1. Introduction

The struggle for computing accurate solution using different grid sizes has increased researchers’ curiosity for developing high order difference schemes. Compact finite-difference scheme is widely used in vast area of computational problems, such as the Helmholtz equations and other elliptic equations [1, 2]. We seek high-accuracy numerical solution of the three-dimensional Helmholtz equation as follows:
(1)
∇
2
u
(
x
,
y
,
z
)
+
l
2
u
(
x
,
y
,
z
)
=
f
(
x
,
y
,
z
)
,
(
x
,
y
,
z
)
∈
Ω
,
where Ω is a cubic solid domain and l is a wave number. The above equation is an elliptic partial differential equation. This equation has broad application in physical phenomena, such as elasticity, electromagnetic waves, potential in time harmonic acoustic and electromagnetic fields, acoustic wave scattering, water wave propagation, noise reduction in silencers, radar scattering, and membrane vibration which are governed in frequency domain. The solution u
(
x
,
y
,
z
) and the forcing function f
(
x
,
y
,
z
) are assumed to have the required continuous partial derivatives up to specific orders and are sufficiently smooth. In this study we analyze finite-difference approximation on uniform grids Δ
x
=
Δ
y
=
Δ
z in x, y, and z directions. In this study we also use a constant value of l to obtain a scheme with sixth order of accuracy [3]. The second-order central-difference operator at x
i, denoted by δ
x
2, is defined as follows:
(2)
δ
x
2
u
i
,
j
,
k
=
u
i
+
1
,
j
,
k
-
2
u
i
,
j
,
k
+
u
i
-
1
,
j
,
k
h
1
2
.
Difference operators δ
y
2 and δ
z
2 are defined similarly. Equation (1) can be discretized at given grid points (
x
i, y
j, and z
k
) as follows:
(3)
δ
x
2
u
i
,
j
,
k
+
δ
y
2
u
i
,
j
,
k
+
δ
z
2
u
i
,
j
,
k
+
l
2
u
i
,
j
,
k
=
f
i
,
j
,
k
+
O
(
h
2
)
,
where O
(
h
2
) denotes the truncated terms on the order O
(
h
1
2
+
h
2
2
+
h
3
2
). Helmholtz equation has been numerically solved by different techniques and different approaches are developed such as the finite-difference method [4], the finite-element method [5], the spectral-element method [6], and the compact finite-difference method [7].In finite-difference methods, the stencil of grid points needs to be enlarged, in order to increase the order of accuracy of approximation, but this is not desirable. Helmholtz equation is extensively solved by finite-element methods, but the disadvantage of finite-element method is the high computational cost; another disadvantage is the pollution effect; that is, its results are less accurate solution at higher wave number for the given nodes per wavelength. In spectral-element method, it shows that it requires fever grid nodes per wavelength as compared to the finite-element methods for the Helmholtz equation [6]. But due to the less sparse resultant matrix, compared to the resulting finite-element matrix, computational time of both methods is the same [6]. Turkel et al. solved Helmholtz equation in 2D and 3D domain for variable wave number l [8]. Many iterative techniques for the Helmholtz equation suffer because of their slow convergence, when high frequencies are required. The struggle for fast iterative methods for high-frequency Helmholtz equations becomes the focus of research.In general to obtain more accurate numerical solution add more nodes and use smaller mesh sizes, which needs more computational time and storage space. Singer and Turkel have conducted an important work in this regard [4]. They used Dirichlet and/or Neumann boundary conditions for the development of a fourth-order compact finite-difference method. Later on a sixth-order finite-difference method was developed by Nabavi et al. for solving Helmholtz equation in one-dimensional and two-dimensional domain with Neumann boundary conditions [7]. It is shown that sixth-order accuracy is the best that can be achieved for Poisson and Helmholtz equation in 2D case [4]. In this work the basic issue discussed is to develop a sixth-order compact finite-difference method for solving three-dimensional Helmholtz equation with Multigrid method. The present study is the first that uses sixth-order compact finite-difference scheme for three dimensions and the designing of specialized Multigrid method. This scheme is different from [1, 7] but the difference is of high order, because the term ∂
6
u
/
∂
x
2
∂
y
2
∂
z
2 in 3 dimensions has no counterpart in 2 dimensions. It is also different when l is kept constant in [8].


## 2. Higher Order Compact Scheme

### 2.1. Two-Dimensional Case

Consider the two-dimensional Helmholtz equation:(4)
u
x
x
(
x
,
y
)
+
u
y
y
(
x
,
y
)
+
l
2
u
=
f
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
.
In order to achieve an appropriate description, from Taylor series expansion, we have
(5)
u
i
+
1
,
j
=
u
i
,
j
+
h
∂
x
u
i
,
j
+
h
2
2
∂
x
2
u
i
,
j
+
h
3
6
∂
x
3
u
i
,
j
+
h
4
24
∂
x
4
u
i
,
j
+
h
5
120
∂
x
5
u
i
,
j
+
h
6
720
∂
x
6
u
i
,
j
+
O
(
h
7
)
,
u
i
-
1
,
j
=
u
i
,
j
-
h
∂
x
u
i
,
j
+
h
2
2
∂
x
2
u
i
,
j
-
h
3
6
∂
x
3
u
i
,
j
+
h
4
24
∂
x
4
u
i
,
j
-
h
5
120
∂
x
5
u
i
,
j
+
h
6
720
∂
x
6
u
i
,
j
+
O
(
h
7
)
,
where the grid space is Δ
x
=
h; adding above expressions (5) and solving for the second-order derivative which gives
(6)
∂
x
2
u
i
,
j
=
u
i
+
1
,
j
-
2
u
i
,
j
+
u
i
-
1
,
j
h
2
-
h
2
12
∂
x
4
u
i
,
j
-
h
4
360
∂
x
6
u
i
,
j
+
O
(
h
6
)
,
using (2) and (6), we have
(7)
δ
x
2
u
i
,
j
=
∂
x
2
u
i
,
j
+
h
2
12
∂
x
4
u
i
,
j
+
h
4
360
∂
x
6
u
i
,
j
+
O
(
h
6
)
;
similarly we can find the approximation for the variable y. Therefore the central difference scheme for Helmholtz equation can be written as follows:
(8)
δ
x
2
u
i
,
j
+
δ
y
2
u
i
,
j
+
l
2
(
u
i
,
j
)
+
α
i
,
j
=
f
i
,
j
+
O
(
h
6
)
,
where
(9)
α
i
,
j
=
-
h
2
12
[
∂
4
u
i
,
j
∂
x
4
+
∂
4
u
i
,
j
∂
y
4
]
-
h
4
360
[
∂
6
u
i
,
j
∂
x
6
+
∂
6
u
i
,
j
∂
y
6
]
+
O
(
h
6
)
.
The higher order derivatives are expressed through mixed derivatives by successively differentiating (4); this process also includes derivative of the forcing function f
(
x
,
y
). Applying the appropriate derivatives of (4), we get
(10)
(
∂
4
u
∂
x
4
)
i
,
j
=
(
∂
2
f
∂
x
2
-
l
2
∂
2
u
∂
x
2
-
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
,
(
∂
4
u
∂
y
4
)
i
,
j
=
(
∂
2
f
∂
y
2
-
l
2
∂
2
u
∂
y
2
-
∂
4
u
∂
y
2
∂
x
2
)
i
,
j
,
and putting the above equation in (9), we get
(11)
α
i
,
j
=
-
h
2
12
[
(
∂
2
f
∂
x
2
+
∂
2
f
∂
y
2
)
i
,
j
-
l
2
(
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
)
i
,
j
-
2
[
∂
4
u
∂
x
2
∂
y
2
]
i
,
j
]
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
]
i
,
j
+
O
(
h
6
)
,
in order to find the fourth-order approximation of [
∂
4
u
/
∂
x
2
∂
y
2
]
i
,
j in above equation, which can be obtained from Taylor series expansion such that
(12)
[
∂
4
u
∂
2
x
∂
2
y
]
i
,
j
=
δ
x
2
δ
y
2
u
i
,
j
-
h
2
12
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
+
O
(
h
4
)
;
putting (12) in (11), we get
(13)
α
i
,
j
=
h
2
12
(
-
∇
2
f
i
,
j
+
2
δ
x
2
δ
y
2
u
i
,
j
+
l
2
f
i
,
j
-
l
4
u
i
,
j
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
5
∂
6
u
∂
x
4
∂
y
2
+
5
∂
6
u
∂
x
2
∂
y
4
]
i
,
j
+
O
(
h
6
)
,
where ∇
2 is laplace operator; for compact sixth-order approximation, differentiating (4), we have
(14)
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
l
2
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
=
∂
4
f
i
,
j
∂
x
2
∂
y
2
,
(15)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
)
i
,
j
=
∇
4
f
i
,
j
-
l
2
(
∂
4
u
∂
x
4
+
∂
4
u
∂
y
4
)
i
,
j
-
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
,
where ∇
4 is biharmonic operator; putting (10) and (14) in (15), we get
(16)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
)
i
,
j
=
∇
4
f
i
,
j
-
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
+
3
l
2
(
∂
4
u
∂
x
2
∂
y
2
-
l
2
∇
2
f
+
l
2
(
-
l
2
u
+
f
)
)
i
,
j
,
and using (14), (16), and (12), we have
(17)
α
i
,
j
=
h
2
12
(
-
∇
2
f
i
,
j
+
2
δ
x
2
δ
y
2
u
i
,
j
+
l
2
f
i
,
j
-
l
4
u
i
,
j
)
-
h
4
360
(
∇
4
f
i
,
j
-
l
2
∇
2
f
i
,
j
+
4
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
-
2
l
2
δ
x
2
δ
y
2
u
i
,
j
-
l
6
u
i
,
j
+
l
4
f
i
,
j
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
)
+
O
(
h
6
)
.
Consequently, the compact sixth-order approximation of the two-dimensional Helmholtz equation can be written in modified form which leads to
(18)
h
2
6
(
1
+
l
2
h
2
30
)
(
δ
x
2
δ
y
2
)
v
i
,
j
+
(
δ
x
2
+
δ
y
2
)
v
i
,
j
+
l
2
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
v
i
,
j
=
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
f
i
,
j
+
(
h
2
12
(
1
-
l
2
h
2
30
)
)
∇
2
f
i
,
j
+
h
4
360
∇
4
f
i
,
j
+
h
4
90
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
,
L
.
H
.
S
=
R
.
H
.
S
,
where v
i
,
j is the approximate solution of (4), such that v
i
,
j
=
u
i
,
j
+
O
(
h
6
). Since
(19)
(
δ
x
2
δ
y
2
)
v
i
,
j
=
1
h
4
[
(
v
i
+
1
,
j
+
1
+
v
i
+
1
,
j
-
1
+
v
i
-
1
,
j
+
1
+
v
i
-
1
,
j
-
1
)
-
2
(
v
i
,
j
+
1
+
v
i
,
j
-
1
+
v
i
+
1
,
j
+
v
i
-
1
,
j
)
+
4
(
v
i
,
j
)
]
,
using (19), (18) can be written as follows:
(20)
L
.
H
.
S
=
a
1
v
i
,
j
+
a
2
L
i
,
j
1
+
a
3
L
i
,
j
2
,
where
(21)
a
1
=
-
10
3
+
l
2
h
2
(
46
45
-
l
2
h
2
12
+
l
4
h
4
360
)
,
a
2
=
2
3
-
l
2
h
2
90
,
a
3
=
1
6
+
l
2
h
2
180
,
L
i
,
j
1
=
v
i
+
1
,
j
+
v
i
-
1
,
j
+
v
i
,
j
+
1
+
v
i
,
j
-
1
,
L
i
,
j
2
=
v
i
+
1
,
j
+
1
+
v
i
+
1
,
j
-
1
+
v
i
-
1
,
j
+
1
+
v
i
-
1
,
j
-
1
.
Also the right-hand side of (18) can be expressed as follows:
(22)
R
.
H
.
S
=
b
11
f
i
,
j
+
b
12
∇
2
f
i
,
j
+
b
13
∇
4
f
i
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
,
=
b
1
f
i
,
j
+
b
2
R
i
,
j
1
+
b
3
R
i
,
j
2
,
where
(23)
b
1
=
h
2
90
(
63
-
13
h
2
l
2
2
+
h
4
l
4
4
)
,
b
2
=
h
2
360
(
23

-

h
2
l
2
)
,
b
3
=
h
2
90
,
R
i
,
j
1
=
f
i
+
1
,
j
+
f
i
-
1
,
j
+
f
i
,
j
+
1
+
f
i
,
j
-
1
,
R
i
,
j
2
=
f
i
+
1
,
j
+
1
+
f
i
+
1
,
j
-
1
+
f
i
-
1
,
j
+
1
+
f
i
-
1
,
j
-
1
;
therefore (18) can be written as follows:
(24)
a
1
v
i
,
j
+
a
2
L
i
,
j
1
+
a
3
L
i
,
j
2
=
b
1
f
i
,
j
+
b
2
R
i
,
j
1
+
b
3
R
i
,
j
2
.
This is the compact sixth-order approximation of two-dimensional Helmholtz equation, which can be written in the form A
U
=
b, where A is a sparse, symmetric, and block tridiagonal matrix.


### 2.2. Three-Dimensional Case

Consider the three-dimensional Helmholtz equation:(25)
u
x
x
(
x
,
y
,
z
)
+
u
y
y
(
x
,
y
,
z
)
+
u
z
z
(
x
,
y
,
z
)
+
l
2
u
(
x
,
y
,
z
)
=
f
(
x
,
y
,
z
)
,
(
x
,
y
,
z
)
∈
Ω
.
In order to achieve an appropriate description, from Taylor series expansion, we have
(26)
u
i
+
1
,
j
,
k
=
u
i
,
j
,
k
+
h
∂
x
u
i
,
j
,
k
+
h
2
2
∂
x
2
u
i
,
j
,
k
+
h
3
6
∂
x
3
u
i
,
j
,
k
+
h
4
24
∂
x
4
u
i
,
j
,
k
+
h
5
120
∂
x
5
u
i
,
j
,
k
+
h
6
720
∂
x
6
u
i
,
j
,
k
+
O
(
h
7
)
,
u
i
-
1
,
j
,
k
=
u
i
,
j
,
k
-
h
∂
x
u
i
,
j
,
k
+
h
2
2
∂
x
2
u
i
,
j
,
k
-
h
3
6
∂
x
3
u
i
,
j
,
k
+
h
4
24
∂
x
4
u
i
,
j
,
k
-
h
5
120
∂
x
5
u
i
,
j
,
k
+
h
6
720
∂
x
6
u
i
,
j
,
k
+
O
(
h
7
)
,
and adding the above expressions and solving for the second derivative which gives
(27)
∂
x
2
u
i
,
j
,
k
=
u
i
+
1
,
j
,
k
-
2
u
i
,
j
,
k
+
u
i
-
1
,
j
,
k
h
2
-
h
2
12
∂
x
4
u
i
,
j
,
k
-
h
4
360
∂
x
6
u
i
,
j
,
k
+
O
(
h
6
)
,
using (2) and (27), we have
(28)
δ
x
2
u
i
,
j
,
k
=
∂
x
2
u
i
,
j
,
k
+
h
2
12
∂
x
4
u
i
,
j
,
k
+
h
4
360
∂
x
6
u
i
,
j
,
k
+
O
(
h
6
)
;
similarly we can find the approximation for the variables y and z. Therefore the central difference scheme for Helmholtz equation can be written as follows:
(29)
δ
x
2
u
i
,
j
,
k
+
δ
y
2
u
i
,
j
,
k
+
δ
z
2
u
i
,
j
,
k
+
l
2
(
u
i
,
j
,
k
)
+
α
i
,
j
,
k
=
f
i
,
j
,
k
+
O
(
h
6
)
,
where
(30)
α
i
,
j
,
k
=
-
h
2
12
[
∂
4
u
∂
x
4
+
∂
4
u
∂
y
4
+
∂
4
u
∂
z
4
]
i
,
j
,
k
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
+
O
(
h
6
)
.
For simplicity we use u
(
x
i
,
y
j
,
z
k
)
=
u
i
,
j
,
k and f
(
x
i
,
y
j
,
z
k
)
=
f
i
,
j
,
k. The higher order derivatives are expressed through mixed derivatives by successively differentiating (25); this process also includes derivative of the forcing function f
(
x
,
y
,
z
). Differentiating (25) twice with respect to x
,
y
,
z and then solving it for ∂
4
u
/
∂
x
4, ∂
4
u
/
∂
y
4, and ∂
4
u
/
∂
z
4, respectively, we have
(31)
(
∂
4
u
∂
x
4
)
i
,
j
,
k
=
(
∂
2
f
∂
x
2
-
l
2
∂
2
u
∂
x
2
-
∂
4
u
∂
x
2
∂
y
2
-
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
,
(
∂
4
u
∂
y
4
)
i
,
j
,
k
=
(
∂
2
f
∂
y
2
-
l
2
∂
2
u
∂
y
2
-
∂
4
u
∂
y
2
∂
x
2
-
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
,
(
∂
4
u
∂
z
4
)
i
,
j
,
k
=
(
∂
2
f
∂
z
2
-
l
2
∂
2
u
∂
z
2
-
∂
4
u
∂
z
2
∂
x
2
-
∂
4
u
∂
z
2
∂
y
2
)
i
,
j
,
k
;
therefore (30) becomes
(32)
α
i
,
j
,
k
=
-
h
2
12
[
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
∇
2
f
i
,
j
,
k
-
l
2
∇
2
u
i
,
j
,
k
-
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
]
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
+
O
(
h
6
)
,
where ∇
2 is laplace operator; for compact sixth-order scheme of 3D problem consider the following discretization:
(33)
[
∂
4
u
∂
x
2
∂
y
2
]
i
,
j
,
k
=
δ
x
2
δ
y
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
,
k
+
O
(
h
4
)
,
[
∂
4
u
∂
x
2
∂
z
2
]
i
,
j
,
k
=
δ
x
2
δ
z
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
z
4
)
i
,
j
,
k
+
O
(
h
4
)
,
[
∂
4
u
∂
y
2
∂
z
2
]
i
,
j
,
k
=
δ
y
2
δ
z
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
y
4
∂
z
2
+
∂
6
u
∂
y
2
∂
z
4
)
i
,
j
,
k
+
O
(
h
4
)
;
now (32) becomes
(34)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
i
,
j
,
k
+
l
2
∇
2
u
i
,
j
,
k
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
i
,
j
,
k
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
-
h
4
72
[
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
z
4
+
∂
6
u
∂
z
4
∂
y
2
+
∂
6
u
∂
z
2
∂
y
4
]
i
,
j
,
k
+
O
(
h
6
)
.Hence(35)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
i
,
j
,
k
+
l
2
∇
2
u
i
,
j
,
k
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
i
,
j
,
k
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
+
5
(
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
z
4
+
∂
6
u
∂
z
2
∂
x
4
+
∂
6
u
∂
z
2
∂
y
4
+
∂
6
u
∂
z
4
∂
y
2
)
]
i
,
j
,
k
+
O
(
h
6
)
.For compact sixth-order approximation, from (25), we have
(36)
∂
4
f
i
,
j
,
k
∂
x
2
∂
y
2
=
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
,
k
∂
4
f
i
,
j
,
k
∂
x
2
∂
z
2
=
(
∂
6
u
∂
z
4
∂
x
2
+
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
∂
4
f
i
,
j
,
k
∂
y
2
∂
z
2
=
(
∂
6
u
∂
y
4
∂
z
2
+
∂
6
u
∂
y
2
∂
z
4
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
z
2
∂
y
2
)
i
,
j
,
k
,
and also
(37)
(
∂
6
u
∂
x
6
)
i
,
j
,
k
=
(
∂
4
f
∂
x
4
-
∂
6
u
∂
x
4
∂
y
2
-
∂
6
u
∂
x
4
∂
z
2
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
x
4
)
i
,
j
,
k
,
(
∂
6
u
∂
y
6
)
i
,
j
,
k
=
(
∂
4
f
∂
y
4
-
∂
6
u
∂
x
2
∂
y
4
-
∂
6
u
∂
z
2
∂
y
4
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
y
4
)
i
,
j
,
k
,
(38)
(
∂
6
u
∂
z
6
)
i
,
j
,
k
=
(
∂
4
f
∂
z
4
-
∂
6
u
i
,
j
,
k
∂
x
2
∂
z
4
-
∂
6
u
∂
y
2
∂
z
4
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
z
4
)
i
,
j
,
k
;
using (31) and (36) in (37), we get
(39)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
)
i
,
j
,
k
=
∇
4
f
i
,
j
,
k
-
l
2
∇
4
u
i
,
j
,
k
-
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
+
l
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
y
2
∂
z
2
+
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
+
3
∂
6
u
i
,
j
,
k
∂
x
2
∂
y
2
∂
z
2
,
and using (19) and (39), we have
(40)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
+
l
2
∇
2
u
+
2
(
∂
2
x
∂
2
y
+
∂
2
x
∂
2
z
+
∂
2
y
∂
2
z
)
u
)
i
,
j
,
k
-
h
4
360
(
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
∇
4
f
-
l
2
∇
4
u
+
4
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
-
12
∂
6
u
∂
x
2
∂
y
2
∂
z
2
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
)
i
,
j
,
k
+
O
(
h
6
)
.
Finally
(41)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
+
l
2
(
f
-
l
2
u
)
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
)
i
,
j
,
k
-
h
4
360
(
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
∇
4
f
-
l
2
∇
4
u
+
4
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
-
12
∂
6
u
∂
x
2
∂
y
2
∂
z
2
-
l
2
∇
2
f
+
l
4
f
-
l
6
u
-
2
l
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
)
i
,
j
,
k
+
O
(
h
6
)
.Consequently, the compact sixth-order approximation of the three-dimensional Helmholtz equation can be written in modified form which leads to the following:(42)
(
δ
x
2
+
δ
y
2
+
δ
z
2
)
v
i
,
j
,
k
+
h
2
6
(
1
+
l
2
h
2
30
)
×
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
v
i
,
j
,
k
+
h
4
30
(
∂
6
v
i
,
j
,
k
∂
x
2
∂
y
2
∂
z
2
)
+
l
2
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
v
i
,
j
,
k
=
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
f
i
,
j
,
k
+
(
h
2
12
(
1
-
l
2
h
2
30
)
)
∇
2
f
i
,
j
,
k
+
h
4
360
∇
4
f
i
,
j
,
k
+
h
4
90
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
.
L
.
H
.
S
=
R
.
H
.
S
,
where v
i
,
j
,
k is the approximate solution to (25) and v
i
,
j
,
k
=
u
i
,
j
,
k
+
O
(
h
6
); also ∇
4 is the biharmonic operator. Since
(43)
(
δ
x
2
+
δ
y
2
+
δ
z
2
)
v
i
,
j
,
k
=
1
h
2
[
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
-
6
(
v
i
,
j
,
k
)
]
,
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
v
i
,
j
,
k
=
1
h
4
[
(
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
-
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
+
1
)
-
4
(
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
)
+
12
(
v
i
,
j
,
k
)
]
,
(
∂
6
∂
x
2
∂
y
2
∂
z
2
)
v
i
,
j
,
k
=
1
h
6
[
(
v
i
+
1
,
j
+
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
-
1
+
v
i
+
1
,
j
+
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
)
-
2
(
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
-
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
+
1
)
-
8
u
i
,
j
,
k
+
4
(
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
)
]
,
using (43) the left-hand side of (42) is
(44)
L
.
H
.
S
=
a
1
v
i
,
j
,
k
+
a
2
L
i
,
j
,
k
1
+
a
3
L
i
,
j
,
k
2
+
a
4
L
i
,
j
,
k
3
,
where
(45)
a
1
=
-
64
15
+
l
2
h
2
(
16
15
-
l
2
h
2
12
+
l
4
h
4
360
)
,
a
2
=
7
15
-
l
2
h
2
45
,
a
3
=
-
1
15
+
h
2
6
+
l
2
h
2
180
,
a
4
=
1
30
,
L
i
,
j
,
k
1
=
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
,
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
,
L
i
,
j
,
k
2
=
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
+
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
-
1
,
L
i
,
j
,
k
3
=
(
v
i
+
1
,
j
+
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
-
1
+
v
i
+
1
,
j
+
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
)
.
Also the right-hand side of (42) can be expressed as follows:
(46)
R
.
H
.
S
=
b
11
f
i
,
j
,
k
+
b
12
∇
2
f
i
,
j
,
k
+
b
13
∇
4
f
i
,
j
,
k
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
=
b
1
f
i
,
j
,
k
+
b
2
R
i
,
j
,
k
1
+
b
3
R
i
,
j
,
k
2
,
where
(47)
b
1
=
h
2
(
111
180
+
h
2
60
-
h
2
l
2
12
+
h
4
l
4
310
)
,
b
2
=
h
2
360
(
47
-
h
2
l
2
)
,
b
3
=
h
2
90
,
R
i
,
j
,
k
1
=
f
i
+
1
,
j
,
k
+
f
i
-
1
,
j
,
k
+
f
i
,
j
+
1
,
k
+
f
i
,
j
-
1
,
k
+
f
i
,
j
,
k
+
1
+
f
i
,
j
,
k
-
1
,
R
i
,
j
,
k
2
=
f
i
+
1
,
j
+
1
,
k
+
f
i
+
1
,
j
-
1
,
k
+
f
i
-
1
,
j
+
1
,
k
+
f
i
-
1
,
j
-
1
,
k
+
f
i
+
1
,
j
,
k
+
1
+
f
i
+
1
,
j
,
k
-
1
+
f
i
-
1
,
j
,
k
+
1
+
f
i
,
j
-
1
,
k
-
1
+
f
i
,
j
+
1
,
k
+
1
+
f
i
,
j
-
1
,
k
+
1
+
f
i
,
j
+
1
,
k
-
1
+
f
i
,
j
-
1
,
k
-
1
.
Thus the compact sixth-order approximation of three-dimensional Helmholtz equation in simplest form is as follows:
(48)
a
1
v
i
,
j
,
k
+
a
2
L
i
,
j
,
k
1
+
a
3
L
i
,
j
,
k
2
+
a
4
L
i
,
j
,
k
3
=
b
1
f
i
,
j
,
k
+
b
2
R
i
,
j
,
k
1
+
b
3
R
i
,
j
,
k
2
.
Equation (48) can be written in the form A
U
=
b. In this equation we assume that the derivative of the right-hand side function f
i
,
j
,
k can be determined numerically. In case where f
i
,
j
,
k is not known analytically, then one can approximate the derivatives with high order finite differences. We need only a fourth-order accurate approximation of ∇
2
f
i
,
j
,
k and a second-order accurate approximation of ∇
4
f
i
,
j
,
k and [
∂
4
f
/
∂
x
2
∂
y
2
+
∂
4
f
/
∂
x
2
∂
z
2
+
∂
4
f
/
∂
y
2
∂
z
2
]
i
,
j
,
k. Fourth-order accurate approximation to ∇
2
f
i
,
j
,
k is obtained by using (
1
-
(
h
2
/
12
)
δ
2
)
-
1
δ
2
f
i
,
j
,
k.


## 2.1. Two-Dimensional Case

Consider the two-dimensional Helmholtz equation:(4)
u
x
x
(
x
,
y
)
+
u
y
y
(
x
,
y
)
+
l
2
u
=
f
(
x
,
y
)
,
(
x
,
y
)
∈
Ω
.
In order to achieve an appropriate description, from Taylor series expansion, we have
(5)
u
i
+
1
,
j
=
u
i
,
j
+
h
∂
x
u
i
,
j
+
h
2
2
∂
x
2
u
i
,
j
+
h
3
6
∂
x
3
u
i
,
j
+
h
4
24
∂
x
4
u
i
,
j
+
h
5
120
∂
x
5
u
i
,
j
+
h
6
720
∂
x
6
u
i
,
j
+
O
(
h
7
)
,
u
i
-
1
,
j
=
u
i
,
j
-
h
∂
x
u
i
,
j
+
h
2
2
∂
x
2
u
i
,
j
-
h
3
6
∂
x
3
u
i
,
j
+
h
4
24
∂
x
4
u
i
,
j
-
h
5
120
∂
x
5
u
i
,
j
+
h
6
720
∂
x
6
u
i
,
j
+
O
(
h
7
)
,
where the grid space is Δ
x
=
h; adding above expressions (5) and solving for the second-order derivative which gives
(6)
∂
x
2
u
i
,
j
=
u
i
+
1
,
j
-
2
u
i
,
j
+
u
i
-
1
,
j
h
2
-
h
2
12
∂
x
4
u
i
,
j
-
h
4
360
∂
x
6
u
i
,
j
+
O
(
h
6
)
,
using (2) and (6), we have
(7)
δ
x
2
u
i
,
j
=
∂
x
2
u
i
,
j
+
h
2
12
∂
x
4
u
i
,
j
+
h
4
360
∂
x
6
u
i
,
j
+
O
(
h
6
)
;
similarly we can find the approximation for the variable y. Therefore the central difference scheme for Helmholtz equation can be written as follows:
(8)
δ
x
2
u
i
,
j
+
δ
y
2
u
i
,
j
+
l
2
(
u
i
,
j
)
+
α
i
,
j
=
f
i
,
j
+
O
(
h
6
)
,
where
(9)
α
i
,
j
=
-
h
2
12
[
∂
4
u
i
,
j
∂
x
4
+
∂
4
u
i
,
j
∂
y
4
]
-
h
4
360
[
∂
6
u
i
,
j
∂
x
6
+
∂
6
u
i
,
j
∂
y
6
]
+
O
(
h
6
)
.
The higher order derivatives are expressed through mixed derivatives by successively differentiating (4); this process also includes derivative of the forcing function f
(
x
,
y
). Applying the appropriate derivatives of (4), we get
(10)
(
∂
4
u
∂
x
4
)
i
,
j
=
(
∂
2
f
∂
x
2
-
l
2
∂
2
u
∂
x
2
-
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
,
(
∂
4
u
∂
y
4
)
i
,
j
=
(
∂
2
f
∂
y
2
-
l
2
∂
2
u
∂
y
2
-
∂
4
u
∂
y
2
∂
x
2
)
i
,
j
,
and putting the above equation in (9), we get
(11)
α
i
,
j
=
-
h
2
12
[
(
∂
2
f
∂
x
2
+
∂
2
f
∂
y
2
)
i
,
j
-
l
2
(
∂
2
u
∂
x
2
+
∂
2
u
∂
y
2
)
i
,
j
-
2
[
∂
4
u
∂
x
2
∂
y
2
]
i
,
j
]
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
]
i
,
j
+
O
(
h
6
)
,
in order to find the fourth-order approximation of [
∂
4
u
/
∂
x
2
∂
y
2
]
i
,
j in above equation, which can be obtained from Taylor series expansion such that
(12)
[
∂
4
u
∂
2
x
∂
2
y
]
i
,
j
=
δ
x
2
δ
y
2
u
i
,
j
-
h
2
12
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
+
O
(
h
4
)
;
putting (12) in (11), we get
(13)
α
i
,
j
=
h
2
12
(
-
∇
2
f
i
,
j
+
2
δ
x
2
δ
y
2
u
i
,
j
+
l
2
f
i
,
j
-
l
4
u
i
,
j
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
5
∂
6
u
∂
x
4
∂
y
2
+
5
∂
6
u
∂
x
2
∂
y
4
]
i
,
j
+
O
(
h
6
)
,
where ∇
2 is laplace operator; for compact sixth-order approximation, differentiating (4), we have
(14)
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
l
2
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
=
∂
4
f
i
,
j
∂
x
2
∂
y
2
,
(15)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
)
i
,
j
=
∇
4
f
i
,
j
-
l
2
(
∂
4
u
∂
x
4
+
∂
4
u
∂
y
4
)
i
,
j
-
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
,
where ∇
4 is biharmonic operator; putting (10) and (14) in (15), we get
(16)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
)
i
,
j
=
∇
4
f
i
,
j
-
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
+
3
l
2
(
∂
4
u
∂
x
2
∂
y
2
-
l
2
∇
2
f
+
l
2
(
-
l
2
u
+
f
)
)
i
,
j
,
and using (14), (16), and (12), we have
(17)
α
i
,
j
=
h
2
12
(
-
∇
2
f
i
,
j
+
2
δ
x
2
δ
y
2
u
i
,
j
+
l
2
f
i
,
j
-
l
4
u
i
,
j
)
-
h
4
360
(
∇
4
f
i
,
j
-
l
2
∇
2
f
i
,
j
+
4
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
-
2
l
2
δ
x
2
δ
y
2
u
i
,
j
-
l
6
u
i
,
j
+
l
4
f
i
,
j
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
)
+
O
(
h
6
)
.
Consequently, the compact sixth-order approximation of the two-dimensional Helmholtz equation can be written in modified form which leads to
(18)
h
2
6
(
1
+
l
2
h
2
30
)
(
δ
x
2
δ
y
2
)
v
i
,
j
+
(
δ
x
2
+
δ
y
2
)
v
i
,
j
+
l
2
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
v
i
,
j
=
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
f
i
,
j
+
(
h
2
12
(
1
-
l
2
h
2
30
)
)
∇
2
f
i
,
j
+
h
4
360
∇
4
f
i
,
j
+
h
4
90
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
,
L
.
H
.
S
=
R
.
H
.
S
,
where v
i
,
j is the approximate solution of (4), such that v
i
,
j
=
u
i
,
j
+
O
(
h
6
). Since
(19)
(
δ
x
2
δ
y
2
)
v
i
,
j
=
1
h
4
[
(
v
i
+
1
,
j
+
1
+
v
i
+
1
,
j
-
1
+
v
i
-
1
,
j
+
1
+
v
i
-
1
,
j
-
1
)
-
2
(
v
i
,
j
+
1
+
v
i
,
j
-
1
+
v
i
+
1
,
j
+
v
i
-
1
,
j
)
+
4
(
v
i
,
j
)
]
,
using (19), (18) can be written as follows:
(20)
L
.
H
.
S
=
a
1
v
i
,
j
+
a
2
L
i
,
j
1
+
a
3
L
i
,
j
2
,
where
(21)
a
1
=
-
10
3
+
l
2
h
2
(
46
45
-
l
2
h
2
12
+
l
4
h
4
360
)
,
a
2
=
2
3
-
l
2
h
2
90
,
a
3
=
1
6
+
l
2
h
2
180
,
L
i
,
j
1
=
v
i
+
1
,
j
+
v
i
-
1
,
j
+
v
i
,
j
+
1
+
v
i
,
j
-
1
,
L
i
,
j
2
=
v
i
+
1
,
j
+
1
+
v
i
+
1
,
j
-
1
+
v
i
-
1
,
j
+
1
+
v
i
-
1
,
j
-
1
.
Also the right-hand side of (18) can be expressed as follows:
(22)
R
.
H
.
S
=
b
11
f
i
,
j
+
b
12
∇
2
f
i
,
j
+
b
13
∇
4
f
i
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
i
,
j
,
=
b
1
f
i
,
j
+
b
2
R
i
,
j
1
+
b
3
R
i
,
j
2
,
where
(23)
b
1
=
h
2
90
(
63
-
13
h
2
l
2
2
+
h
4
l
4
4
)
,
b
2
=
h
2
360
(
23

-

h
2
l
2
)
,
b
3
=
h
2
90
,
R
i
,
j
1
=
f
i
+
1
,
j
+
f
i
-
1
,
j
+
f
i
,
j
+
1
+
f
i
,
j
-
1
,
R
i
,
j
2
=
f
i
+
1
,
j
+
1
+
f
i
+
1
,
j
-
1
+
f
i
-
1
,
j
+
1
+
f
i
-
1
,
j
-
1
;
therefore (18) can be written as follows:
(24)
a
1
v
i
,
j
+
a
2
L
i
,
j
1
+
a
3
L
i
,
j
2
=
b
1
f
i
,
j
+
b
2
R
i
,
j
1
+
b
3
R
i
,
j
2
.
This is the compact sixth-order approximation of two-dimensional Helmholtz equation, which can be written in the form A
U
=
b, where A is a sparse, symmetric, and block tridiagonal matrix.


## 2.2. Three-Dimensional Case

Consider the three-dimensional Helmholtz equation:(25)
u
x
x
(
x
,
y
,
z
)
+
u
y
y
(
x
,
y
,
z
)
+
u
z
z
(
x
,
y
,
z
)
+
l
2
u
(
x
,
y
,
z
)
=
f
(
x
,
y
,
z
)
,
(
x
,
y
,
z
)
∈
Ω
.
In order to achieve an appropriate description, from Taylor series expansion, we have
(26)
u
i
+
1
,
j
,
k
=
u
i
,
j
,
k
+
h
∂
x
u
i
,
j
,
k
+
h
2
2
∂
x
2
u
i
,
j
,
k
+
h
3
6
∂
x
3
u
i
,
j
,
k
+
h
4
24
∂
x
4
u
i
,
j
,
k
+
h
5
120
∂
x
5
u
i
,
j
,
k
+
h
6
720
∂
x
6
u
i
,
j
,
k
+
O
(
h
7
)
,
u
i
-
1
,
j
,
k
=
u
i
,
j
,
k
-
h
∂
x
u
i
,
j
,
k
+
h
2
2
∂
x
2
u
i
,
j
,
k
-
h
3
6
∂
x
3
u
i
,
j
,
k
+
h
4
24
∂
x
4
u
i
,
j
,
k
-
h
5
120
∂
x
5
u
i
,
j
,
k
+
h
6
720
∂
x
6
u
i
,
j
,
k
+
O
(
h
7
)
,
and adding the above expressions and solving for the second derivative which gives
(27)
∂
x
2
u
i
,
j
,
k
=
u
i
+
1
,
j
,
k
-
2
u
i
,
j
,
k
+
u
i
-
1
,
j
,
k
h
2
-
h
2
12
∂
x
4
u
i
,
j
,
k
-
h
4
360
∂
x
6
u
i
,
j
,
k
+
O
(
h
6
)
,
using (2) and (27), we have
(28)
δ
x
2
u
i
,
j
,
k
=
∂
x
2
u
i
,
j
,
k
+
h
2
12
∂
x
4
u
i
,
j
,
k
+
h
4
360
∂
x
6
u
i
,
j
,
k
+
O
(
h
6
)
;
similarly we can find the approximation for the variables y and z. Therefore the central difference scheme for Helmholtz equation can be written as follows:
(29)
δ
x
2
u
i
,
j
,
k
+
δ
y
2
u
i
,
j
,
k
+
δ
z
2
u
i
,
j
,
k
+
l
2
(
u
i
,
j
,
k
)
+
α
i
,
j
,
k
=
f
i
,
j
,
k
+
O
(
h
6
)
,
where
(30)
α
i
,
j
,
k
=
-
h
2
12
[
∂
4
u
∂
x
4
+
∂
4
u
∂
y
4
+
∂
4
u
∂
z
4
]
i
,
j
,
k
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
+
O
(
h
6
)
.
For simplicity we use u
(
x
i
,
y
j
,
z
k
)
=
u
i
,
j
,
k and f
(
x
i
,
y
j
,
z
k
)
=
f
i
,
j
,
k. The higher order derivatives are expressed through mixed derivatives by successively differentiating (25); this process also includes derivative of the forcing function f
(
x
,
y
,
z
). Differentiating (25) twice with respect to x
,
y
,
z and then solving it for ∂
4
u
/
∂
x
4, ∂
4
u
/
∂
y
4, and ∂
4
u
/
∂
z
4, respectively, we have
(31)
(
∂
4
u
∂
x
4
)
i
,
j
,
k
=
(
∂
2
f
∂
x
2
-
l
2
∂
2
u
∂
x
2
-
∂
4
u
∂
x
2
∂
y
2
-
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
,
(
∂
4
u
∂
y
4
)
i
,
j
,
k
=
(
∂
2
f
∂
y
2
-
l
2
∂
2
u
∂
y
2
-
∂
4
u
∂
y
2
∂
x
2
-
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
,
(
∂
4
u
∂
z
4
)
i
,
j
,
k
=
(
∂
2
f
∂
z
2
-
l
2
∂
2
u
∂
z
2
-
∂
4
u
∂
z
2
∂
x
2
-
∂
4
u
∂
z
2
∂
y
2
)
i
,
j
,
k
;
therefore (30) becomes
(32)
α
i
,
j
,
k
=
-
h
2
12
[
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
∇
2
f
i
,
j
,
k
-
l
2
∇
2
u
i
,
j
,
k
-
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
i
,
j
,
k
]
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
+
O
(
h
6
)
,
where ∇
2 is laplace operator; for compact sixth-order scheme of 3D problem consider the following discretization:
(33)
[
∂
4
u
∂
x
2
∂
y
2
]
i
,
j
,
k
=
δ
x
2
δ
y
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
)
i
,
j
,
k
+
O
(
h
4
)
,
[
∂
4
u
∂
x
2
∂
z
2
]
i
,
j
,
k
=
δ
x
2
δ
z
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
z
4
)
i
,
j
,
k
+
O
(
h
4
)
,
[
∂
4
u
∂
y
2
∂
z
2
]
i
,
j
,
k
=
δ
y
2
δ
z
2
u
i
,
j
,
k
-
h
2
12
(
∂
6
u
∂
y
4
∂
z
2
+
∂
6
u
∂
y
2
∂
z
4
)
i
,
j
,
k
+
O
(
h
4
)
;
now (32) becomes
(34)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
i
,
j
,
k
+
l
2
∇
2
u
i
,
j
,
k
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
i
,
j
,
k
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
]
i
,
j
,
k
-
h
4
72
[
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
z
4
+
∂
6
u
∂
z
4
∂
y
2
+
∂
6
u
∂
z
2
∂
y
4
]
i
,
j
,
k
+
O
(
h
6
)
.Hence(35)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
i
,
j
,
k
+
l
2
∇
2
u
i
,
j
,
k
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
i
,
j
,
k
)
-
h
4
360
[
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
+
5
(
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
z
4
+
∂
6
u
∂
z
2
∂
x
4
+
∂
6
u
∂
z
2
∂
y
4
+
∂
6
u
∂
z
4
∂
y
2
)
]
i
,
j
,
k
+
O
(
h
6
)
.For compact sixth-order approximation, from (25), we have
(36)
∂
4
f
i
,
j
,
k
∂
x
2
∂
y
2
=
(
∂
6
u
∂
x
4
∂
y
2
+
∂
6
u
∂
x
2
∂
y
4
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
x
2
∂
y
2
)
i
,
j
,
k
∂
4
f
i
,
j
,
k
∂
x
2
∂
z
2
=
(
∂
6
u
∂
z
4
∂
x
2
+
∂
6
u
∂
x
4
∂
z
2
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
∂
4
f
i
,
j
,
k
∂
y
2
∂
z
2
=
(
∂
6
u
∂
y
4
∂
z
2
+
∂
6
u
∂
y
2
∂
z
4
+
∂
6
u
∂
x
2
∂
y
2
∂
z
2
+
l
2
∂
4
u
∂
z
2
∂
y
2
)
i
,
j
,
k
,
and also
(37)
(
∂
6
u
∂
x
6
)
i
,
j
,
k
=
(
∂
4
f
∂
x
4
-
∂
6
u
∂
x
4
∂
y
2
-
∂
6
u
∂
x
4
∂
z
2
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
x
4
)
i
,
j
,
k
,
(
∂
6
u
∂
y
6
)
i
,
j
,
k
=
(
∂
4
f
∂
y
4
-
∂
6
u
∂
x
2
∂
y
4
-
∂
6
u
∂
z
2
∂
y
4
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
y
4
)
i
,
j
,
k
,
(38)
(
∂
6
u
∂
z
6
)
i
,
j
,
k
=
(
∂
4
f
∂
z
4
-
∂
6
u
i
,
j
,
k
∂
x
2
∂
z
4
-
∂
6
u
∂
y
2
∂
z
4
)
i
,
j
,
k
hhhhhhhhhhh
-
l
2
(
∂
4
u
∂
z
4
)
i
,
j
,
k
;
using (31) and (36) in (37), we get
(39)
(
∂
6
u
∂
x
6
+
∂
6
u
∂
y
6
+
∂
6
u
∂
z
6
)
i
,
j
,
k
=
∇
4
f
i
,
j
,
k
-
l
2
∇
4
u
i
,
j
,
k
-
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
+
l
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
y
2
∂
z
2
+
∂
4
u
∂
x
2
∂
z
2
)
i
,
j
,
k
+
3
∂
6
u
i
,
j
,
k
∂
x
2
∂
y
2
∂
z
2
,
and using (19) and (39), we have
(40)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
+
l
2
∇
2
u
+
2
(
∂
2
x
∂
2
y
+
∂
2
x
∂
2
z
+
∂
2
y
∂
2
z
)
u
)
i
,
j
,
k
-
h
4
360
(
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
∇
4
f
-
l
2
∇
4
u
+
4
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
-
12
∂
6
u
∂
x
2
∂
y
2
∂
z
2
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
)
i
,
j
,
k
+
O
(
h
6
)
.
Finally
(41)
α
i
,
j
,
k
=
h
2
12
(
-
∇
2
f
+
l
2
(
f
-
l
2
u
)
+
2
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
u
)
i
,
j
,
k
-
h
4
360
(
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
∇
4
f
-
l
2
∇
4
u
+
4
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
-
12
∂
6
u
∂
x
2
∂
y
2
∂
z
2
-
l
2
∇
2
f
+
l
4
f
-
l
6
u
-
2
l
2
(
∂
4
u
∂
x
2
∂
y
2
+
∂
4
u
∂
x
2
∂
z
2
+
∂
4
u
∂
y
2
∂
z
2
)
)
i
,
j
,
k
+
O
(
h
6
)
.Consequently, the compact sixth-order approximation of the three-dimensional Helmholtz equation can be written in modified form which leads to the following:(42)
(
δ
x
2
+
δ
y
2
+
δ
z
2
)
v
i
,
j
,
k
+
h
2
6
(
1
+
l
2
h
2
30
)
×
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
v
i
,
j
,
k
+
h
4
30
(
∂
6
v
i
,
j
,
k
∂
x
2
∂
y
2
∂
z
2
)
+
l
2
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
v
i
,
j
,
k
=
(
1
-
l
2
h
2
12
+
l
4
h
4
360
)
f
i
,
j
,
k
+
(
h
2
12
(
1
-
l
2
h
2
30
)
)
∇
2
f
i
,
j
,
k
+
h
4
360
∇
4
f
i
,
j
,
k
+
h
4
90
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
.
L
.
H
.
S
=
R
.
H
.
S
,
where v
i
,
j
,
k is the approximate solution to (25) and v
i
,
j
,
k
=
u
i
,
j
,
k
+
O
(
h
6
); also ∇
4 is the biharmonic operator. Since
(43)
(
δ
x
2
+
δ
y
2
+
δ
z
2
)
v
i
,
j
,
k
=
1
h
2
[
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
-
6
(
v
i
,
j
,
k
)
]
,
(
δ
x
2
δ
y
2
+
δ
x
2
δ
z
2
+
δ
y
2
δ
z
2
)
v
i
,
j
,
k
=
1
h
4
[
(
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
-
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
+
1
)
-
4
(
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
)
+
12
(
v
i
,
j
,
k
)
]
,
(
∂
6
∂
x
2
∂
y
2
∂
z
2
)
v
i
,
j
,
k
=
1
h
6
[
(
v
i
+
1
,
j
+
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
-
1
+
v
i
+
1
,
j
+
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
)
-
2
(
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
-
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
+
1
)
-
8
u
i
,
j
,
k
+
4
(
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
)
]
,
using (43) the left-hand side of (42) is
(44)
L
.
H
.
S
=
a
1
v
i
,
j
,
k
+
a
2
L
i
,
j
,
k
1
+
a
3
L
i
,
j
,
k
2
+
a
4
L
i
,
j
,
k
3
,
where
(45)
a
1
=
-
64
15
+
l
2
h
2
(
16
15
-
l
2
h
2
12
+
l
4
h
4
360
)
,
a
2
=
7
15
-
l
2
h
2
45
,
a
3
=
-
1
15
+
h
2
6
+
l
2
h
2
180
,
a
4
=
1
30
,
L
i
,
j
,
k
1
=
v
i
+
1
,
j
,
k
+
v
i
-
1
,
j
,
k
+
v
i
,
j
+
1
,
k
,
+
v
i
,
j
-
1
,
k
+
v
i
,
j
,
k
+
1
+
v
i
,
j
,
k
-
1
,
L
i
,
j
,
k
2
=
v
i
+
1
,
j
+
1
,
k
+
v
i
+
1
,
j
-
1
,
k
+
v
i
-
1
,
j
+
1
,
k
+
v
i
-
1
,
j
-
1
,
k
+
v
i
+
1
,
j
,
k
+
1
+
v
i
+
1
,
j
,
k
-
1
+
v
i
-
1
,
j
,
k
+
1
+
v
i
-
1
,
j
,
k
-
1
+
v
i
,
j
+
1
,
k
+
1
+
v
i
,
j
-
1
,
k
+
1
+
v
i
,
j
+
1
,
k
-
1
+
v
i
,
j
-
1
,
k
-
1
,
L
i
,
j
,
k
3
=
(
v
i
+
1
,
j
+
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
+
1
+
v
i
+
1
,
j
-
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
+
1
+
v
i
-
1
,
j
-
1
,
k
-
1
+
v
i
+
1
,
j
+
1
,
k
-
1
+
v
i
-
1
,
j
+
1
,
k
+
1
)
.
Also the right-hand side of (42) can be expressed as follows:
(46)
R
.
H
.
S
=
b
11
f
i
,
j
,
k
+
b
12
∇
2
f
i
,
j
,
k
+
b
13
∇
4
f
i
,
j
,
k
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
+
∂
4
f
∂
x
2
∂
z
2
+
∂
4
f
∂
y
2
∂
z
2
]
i
,
j
,
k
=
b
1
f
i
,
j
,
k
+
b
2
R
i
,
j
,
k
1
+
b
3
R
i
,
j
,
k
2
,
where
(47)
b
1
=
h
2
(
111
180
+
h
2
60
-
h
2
l
2
12
+
h
4
l
4
310
)
,
b
2
=
h
2
360
(
47
-
h
2
l
2
)
,
b
3
=
h
2
90
,
R
i
,
j
,
k
1
=
f
i
+
1
,
j
,
k
+
f
i
-
1
,
j
,
k
+
f
i
,
j
+
1
,
k
+
f
i
,
j
-
1
,
k
+
f
i
,
j
,
k
+
1
+
f
i
,
j
,
k
-
1
,
R
i
,
j
,
k
2
=
f
i
+
1
,
j
+
1
,
k
+
f
i
+
1
,
j
-
1
,
k
+
f
i
-
1
,
j
+
1
,
k
+
f
i
-
1
,
j
-
1
,
k
+
f
i
+
1
,
j
,
k
+
1
+
f
i
+
1
,
j
,
k
-
1
+
f
i
-
1
,
j
,
k
+
1
+
f
i
,
j
-
1
,
k
-
1
+
f
i
,
j
+
1
,
k
+
1
+
f
i
,
j
-
1
,
k
+
1
+
f
i
,
j
+
1
,
k
-
1
+
f
i
,
j
-
1
,
k
-
1
.
Thus the compact sixth-order approximation of three-dimensional Helmholtz equation in simplest form is as follows:
(48)
a
1
v
i
,
j
,
k
+
a
2
L
i
,
j
,
k
1
+
a
3
L
i
,
j
,
k
2
+
a
4
L
i
,
j
,
k
3
=
b
1
f
i
,
j
,
k
+
b
2
R
i
,
j
,
k
1
+
b
3
R
i
,
j
,
k
2
.
Equation (48) can be written in the form A
U
=
b. In this equation we assume that the derivative of the right-hand side function f
i
,
j
,
k can be determined numerically. In case where f
i
,
j
,
k is not known analytically, then one can approximate the derivatives with high order finite differences. We need only a fourth-order accurate approximation of ∇
2
f
i
,
j
,
k and a second-order accurate approximation of ∇
4
f
i
,
j
,
k and [
∂
4
f
/
∂
x
2
∂
y
2
+
∂
4
f
/
∂
x
2
∂
z
2
+
∂
4
f
/
∂
y
2
∂
z
2
]
i
,
j
,
k. Fourth-order accurate approximation to ∇
2
f
i
,
j
,
k is obtained by using (
1
-
(
h
2
/
12
)
δ
2
)
-
1
δ
2
f
i
,
j
,
k.


## 3. Boundary Condition

When a Dirichlet boundary condition is applied, then the formula in (18) can be used for all interior points. In case of Neumann boundary condition, we develop a sixth-order accurate method for two- and three-dimensional cases; that is, u
x
=
g
(
y
) and u
x
=
g
(
y
,
z
). Similar formulas hold in the other directions. To be specific, we consider the coordinate line i
=
0, in (18) and introduce a ghost point i
=
-
1. At the boundary i
=
0, we specify both the Helmholtz equation and the Neumann boundary condition.


### 3.1. Two-Dimensional Case

Considering (4), using (5), we have
(49)
δ
x
0
u
i
=
∂
x
u
i
,
j
+
h
2
6
∂
x
3
u
i
,
j
+
h
4
120
∂
x
5
u
i
,
j
+
O
(
h
6
)
.
From (4), we have
(50)
(
∂
3
u
∂
x
3
)
i
,
j
=
[
(
∂
f
∂
x
-
∂
3
u
∂
x
∂
y
2
)
-
l
2
(
∂
u
∂
x
)
]
i
,
j
,
(
∂
5
u
∂
x
5
)
i
,
j
=
[
(
∂
3
f
∂
x
3
-
∂
5
u
∂
x
3
∂
y
2
)
-
l
2
(
∂
3
u
∂
x
3
)
]
i
,
j
.
In above equation (50), we need fourth-order approximation of (
∂
3
u
/
∂
x
∂
y
2
); that is,
(51)
(
∂
3
u
∂
x
∂
y
2
)
i
,
j
=
δ
x
δ
y
2
u
i
,
j
-
h
2
6
(
∂
5
u
∂
x
∂
y
4
+
2
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
+
O
(
h
4
)
.
Making use of (51) in (50), we have
(52)
(
∂
3
u
∂
x
3
)
i
,
j
=
h
2
12
(
∂
5
u
∂
x
∂
y
4
+
2
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
-
δ
x
δ
y
2
u
i
,
j
-
l
2
(
∂
u
∂
x
)
i
,
j
+
(
∂
3
f
∂
x
3
)
i
,
j
+
O
(
h
4
)
.
Now the second-order approximation of (
∂
3
u
/
∂
x
3
) is
(53)
(
∂
3
u
∂
x
3
)
i
,
j
=
(
∂
f
∂
x
)
i
,
j
-
δ
x
δ
y
2
u
i
,
j
-
l
2
(
∂
u
∂
x
)
i
,
j
+
O
(
h
2
)
,
and also the second-order approximation of (
∂
5
u
/
∂
x
5
) is
(54)
(
∂
5
u
∂
x
5
)
i
,
j
=
(
∂
3
f
∂
x
3
)
i
,
j
-
l
2
(
∂
f
∂
x
)
i
,
j
+
l
4
(
∂
u
∂
x
)
i
,
j
-
l
2
δ
x
δ
y
2
u
i
,
j
-
(
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
+
O
(
h
2
)
.
Using the derivatives of (4), we have
(55)
(
∂
5
u
∂
x
3
∂
y
2
+
∂
5
u
∂
x
∂
y
4
)
i
,
j
=
(
∂
3
f
∂
x
3
-
l
2
(
∂
3
u
∂
x
∂
y
2
)
)
i
,
j
.
Using (52), (54), and (55) in (49), we have
(56)
δ
x
0
u
i
-
l
4
h
4
120
(
∂
u
∂
x
)
i
,
j
+
h
2
6
(
1
+
l
2
h
2
30
)
δ
x
δ
y
2
u
i
,
j
=
(
1
-
l
2
h
2
6
(
∂
u
∂
x
)
i
,
j
)
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
i
,
j
+
h
4
120
(
∂
3
u
∂
x
3
)
i
,
j
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
i
,
j
+
O
(
h
6
)
,
for (
∂
u
/
∂
x
); in (56) the approximation used is
(57)
(
∂
u
∂
x
)
i
,
j
=
δ
x
0
u
i
+
β
h
2
δ
x
δ
y
2
u
i
,
j
,
where β is arbitrary constant; using (57) in (56), we have
(58)
L
1
*
(
v
1
,
j
-
v
-
1
,
j
)
+
L
2
*
(
v
1
,
j
+
1
+
v
1
,
j
-
1
+
v
-
1
,
j
+
1
+
v
-
1
,
j
-
1
)
=
[
(
1
-
l
2
h
2
6
)
g
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
+
h
4
120
(
∂
3
u
∂
x
3
)
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
]
i
,
j
,
where g
=
g
(
y
), L
1
*
=
1
-
(
l
4
h
4
/
120
)
(
1
-
2
β
)  −  (
1
/
6
h
)
(
1
+
l
2
h
2
/
30
), and L
2
*
=
-
(
l
4
h
4
/
120
)
(
-
β
)
-
(
1
/
12
h
)
(
1
+
l
2
h
2
/
30
). Putting i
=
0 in (24), we have
(59)
a
1
v
0
,
j
+
a
2
(
v
1
,
j
+
v
-
1
,
j
+
v
0
,
j
+
1
+
v
0
,
j
-
1
)
+
a
3
(
v
1
,
j
+
1
+
v
1
,
j
-
1
+
v
-
1
,
j
+
1
+
v
-
1
,
j
-
1
)
=
b
11
f
0
,
j
+
b
12
∇
2
f
0
,
j
+
b
13
∇
4
f
0
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
0
,
j
.
Multiplying (58) by β and adding (59), we get formula for boundary nodes
(60)
a
1
v
0
,
j
+
a
2
(
2
v
1
,
j
+
v
0
,
j
+
1
+
v
0
,
j
-
1
)
+
a
3
(
v
1
,
j
+
1
+
v
1
,
j
-
1
)
=
b
11
f
0
,
j
+
b
12
∇
2
f
0
,
j
+
b
13
∇
4
f
0
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
0
,
j
+
β
[
(
1
-
l
2
h
2
6
)
g
y
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
0
,
j
+
h
4
120
(
∂
3
u
∂
x
3
)
0
,
j
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
0
,
j
]
,
where β
=
a
2
/
L
1
*
=
a
3
/
L
2
*
=
120
/
(
120
-
l
4
h
4
) in order to eliminate all the elements with i
=
-
1.


### 3.2. Three-Dimensional Case

Consider (25), by using (26), we have
(61)
u
i
+
1
-
u
i
-
1
2
h
=
∂
x
u
+
h
2
6
∂
x
3
u
+
h
4
120
∂
x
5
u
+
O
(
h
6
)
.
To be specific, we consider the coordinate line i
=
0 and introduce a ghost point i
=
-
1; from (25), we have
(62)
(
∂
3
u
∂
x
3
)
=
-
(
∂
3
u
∂
x
∂
y
2
)
-
(
∂
3
u
∂
x
∂
z
2
)
-
l
2
(
∂
u
∂
x
)
+
∂
f
∂
x
,
(63)
(
∂
5
u
∂
x
5
)
=
-
(
∂
5
u
∂
x
3
∂
y
2
)
-
(
∂
5
u
∂
x
3
∂
z
2
)
-
l
2
(
∂
3
u
∂
x
3
)
+
∂
3
f
∂
x
3
.
Using (62) in (63), we have
(64)
(
∂
5
u
∂
x
5
)
=
(
∂
5
u
∂
x
3
∂
y
2
+
∂
5
u
∂
x
3
∂
z
2
)
-
l
2
(
-
∂
3
u
∂
x
∂
y
2
-
∂
3
u
∂
x
∂
z
2
-
l
2
(
∂
u
∂
x
)
+
∂
f
∂
x
)
+
∂
3
f
∂
x
3
=
(
∂
5
u
∂
x
∂
y
4
+
∂
5
u
∂
x
∂
z
4
)
+
2
∂
5
u
∂
x
∂
y
2
∂
z
2
+
2
l
2
(
∂
3
u
∂
x
∂
y
2
+
∂
3
u
∂
x
∂
z
2
)
+
l
4
(
∂
u
∂
x
)
-
l
2
∂
f
∂
x
+
∂
3
f
∂
x
3
-
∂
3
f
∂
x
∂
y
2
-
∂
3
f
∂
x
∂
z
2
.
Now using center difference discretization at (
0
,
j
,
k
), (25) will become
(65)
g
-
h
2
6
[
∂
2
g
∂
y
2
+
∂
2
g
∂
z
2
+
l
2
(
g
)
-
∂
f
∂
x
]
+
h
4
120
×
[
∂
4
g
∂
y
4
+
∂
4
g
∂
z
4
+
2
∂
4
g
∂
y
2
∂
z
2
+
2
l
2
(
∂
2
g
∂
y
2
+
∂
2
g
∂
z
2
)
+
l
4
(
g
)
-
l
2
∂
f
∂
x
+
∂
3
f
∂
x
3
-
∂
3
f
∂
x
∂
y
2
-
∂
3
f
∂
x
∂
z
2
]
.
Equations (60) and (65) yield that for Neumann boundary points the stencil is of five points for two-dimensional case while in three-dimensional case the stencil has seven points. The matrix is inverted including the extra line i
=
-
1. Applying Neumann boundary conditions accuracy and CPU timing remain the same in all examples.Different methods are used to solve (48) such as finite-difference scheme, fourth-order compact finite-difference scheme, and LU decomposition. We are solving this equation by Multigrid method with sixth-order compact finite-difference scheme using Gauss-Seidel method as a smoother.


## 3.1. Two-Dimensional Case

Considering (4), using (5), we have
(49)
δ
x
0
u
i
=
∂
x
u
i
,
j
+
h
2
6
∂
x
3
u
i
,
j
+
h
4
120
∂
x
5
u
i
,
j
+
O
(
h
6
)
.
From (4), we have
(50)
(
∂
3
u
∂
x
3
)
i
,
j
=
[
(
∂
f
∂
x
-
∂
3
u
∂
x
∂
y
2
)
-
l
2
(
∂
u
∂
x
)
]
i
,
j
,
(
∂
5
u
∂
x
5
)
i
,
j
=
[
(
∂
3
f
∂
x
3
-
∂
5
u
∂
x
3
∂
y
2
)
-
l
2
(
∂
3
u
∂
x
3
)
]
i
,
j
.
In above equation (50), we need fourth-order approximation of (
∂
3
u
/
∂
x
∂
y
2
); that is,
(51)
(
∂
3
u
∂
x
∂
y
2
)
i
,
j
=
δ
x
δ
y
2
u
i
,
j
-
h
2
6
(
∂
5
u
∂
x
∂
y
4
+
2
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
+
O
(
h
4
)
.
Making use of (51) in (50), we have
(52)
(
∂
3
u
∂
x
3
)
i
,
j
=
h
2
12
(
∂
5
u
∂
x
∂
y
4
+
2
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
-
δ
x
δ
y
2
u
i
,
j
-
l
2
(
∂
u
∂
x
)
i
,
j
+
(
∂
3
f
∂
x
3
)
i
,
j
+
O
(
h
4
)
.
Now the second-order approximation of (
∂
3
u
/
∂
x
3
) is
(53)
(
∂
3
u
∂
x
3
)
i
,
j
=
(
∂
f
∂
x
)
i
,
j
-
δ
x
δ
y
2
u
i
,
j
-
l
2
(
∂
u
∂
x
)
i
,
j
+
O
(
h
2
)
,
and also the second-order approximation of (
∂
5
u
/
∂
x
5
) is
(54)
(
∂
5
u
∂
x
5
)
i
,
j
=
(
∂
3
f
∂
x
3
)
i
,
j
-
l
2
(
∂
f
∂
x
)
i
,
j
+
l
4
(
∂
u
∂
x
)
i
,
j
-
l
2
δ
x
δ
y
2
u
i
,
j
-
(
∂
5
u
∂
x
3
∂
y
2
)
i
,
j
+
O
(
h
2
)
.
Using the derivatives of (4), we have
(55)
(
∂
5
u
∂
x
3
∂
y
2
+
∂
5
u
∂
x
∂
y
4
)
i
,
j
=
(
∂
3
f
∂
x
3
-
l
2
(
∂
3
u
∂
x
∂
y
2
)
)
i
,
j
.
Using (52), (54), and (55) in (49), we have
(56)
δ
x
0
u
i
-
l
4
h
4
120
(
∂
u
∂
x
)
i
,
j
+
h
2
6
(
1
+
l
2
h
2
30
)
δ
x
δ
y
2
u
i
,
j
=
(
1
-
l
2
h
2
6
(
∂
u
∂
x
)
i
,
j
)
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
i
,
j
+
h
4
120
(
∂
3
u
∂
x
3
)
i
,
j
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
i
,
j
+
O
(
h
6
)
,
for (
∂
u
/
∂
x
); in (56) the approximation used is
(57)
(
∂
u
∂
x
)
i
,
j
=
δ
x
0
u
i
+
β
h
2
δ
x
δ
y
2
u
i
,
j
,
where β is arbitrary constant; using (57) in (56), we have
(58)
L
1
*
(
v
1
,
j
-
v
-
1
,
j
)
+
L
2
*
(
v
1
,
j
+
1
+
v
1
,
j
-
1
+
v
-
1
,
j
+
1
+
v
-
1
,
j
-
1
)
=
[
(
1
-
l
2
h
2
6
)
g
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
+
h
4
120
(
∂
3
u
∂
x
3
)
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
]
i
,
j
,
where g
=
g
(
y
), L
1
*
=
1
-
(
l
4
h
4
/
120
)
(
1
-
2
β
)  −  (
1
/
6
h
)
(
1
+
l
2
h
2
/
30
), and L
2
*
=
-
(
l
4
h
4
/
120
)
(
-
β
)
-
(
1
/
12
h
)
(
1
+
l
2
h
2
/
30
). Putting i
=
0 in (24), we have
(59)
a
1
v
0
,
j
+
a
2
(
v
1
,
j
+
v
-
1
,
j
+
v
0
,
j
+
1
+
v
0
,
j
-
1
)
+
a
3
(
v
1
,
j
+
1
+
v
1
,
j
-
1
+
v
-
1
,
j
+
1
+
v
-
1
,
j
-
1
)
=
b
11
f
0
,
j
+
b
12
∇
2
f
0
,
j
+
b
13
∇
4
f
0
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
0
,
j
.
Multiplying (58) by β and adding (59), we get formula for boundary nodes
(60)
a
1
v
0
,
j
+
a
2
(
2
v
1
,
j
+
v
0
,
j
+
1
+
v
0
,
j
-
1
)
+
a
3
(
v
1
,
j
+
1
+
v
1
,
j
-
1
)
=
b
11
f
0
,
j
+
b
12
∇
2
f
0
,
j
+
b
13
∇
4
f
0
,
j
+
b
14
[
∂
4
f
∂
x
2
∂
y
2
]
0
,
j
+
β
[
(
1
-
l
2
h
2
6
)
g
y
+
h
2
6
(
1
-
l
2
h
2
20
)
(
∂
f
∂
x
)
0
,
j
+
h
4
120
(
∂
3
u
∂
x
3
)
0
,
j
+
h
4
72
(
∂
3
f
∂
x
∂
y
2
)
0
,
j
]
,
where β
=
a
2
/
L
1
*
=
a
3
/
L
2
*
=
120
/
(
120
-
l
4
h
4
) in order to eliminate all the elements with i
=
-
1.


## 3.2. Three-Dimensional Case

Consider (25), by using (26), we have
(61)
u
i
+
1
-
u
i
-
1
2
h
=
∂
x
u
+
h
2
6
∂
x
3
u
+
h
4
120
∂
x
5
u
+
O
(
h
6
)
.
To be specific, we consider the coordinate line i
=
0 and introduce a ghost point i
=
-
1; from (25), we have
(62)
(
∂
3
u
∂
x
3
)
=
-
(
∂
3
u
∂
x
∂
y
2
)
-
(
∂
3
u
∂
x
∂
z
2
)
-
l
2
(
∂
u
∂
x
)
+
∂
f
∂
x
,
(63)
(
∂
5
u
∂
x
5
)
=
-
(
∂
5
u
∂
x
3
∂
y
2
)
-
(
∂
5
u
∂
x
3
∂
z
2
)
-
l
2
(
∂
3
u
∂
x
3
)
+
∂
3
f
∂
x
3
.
Using (62) in (63), we have
(64)
(
∂
5
u
∂
x
5
)
=
(
∂
5
u
∂
x
3
∂
y
2
+
∂
5
u
∂
x
3
∂
z
2
)
-
l
2
(
-
∂
3
u
∂
x
∂
y
2
-
∂
3
u
∂
x
∂
z
2
-
l
2
(
∂
u
∂
x
)
+
∂
f
∂
x
)
+
∂
3
f
∂
x
3
=
(
∂
5
u
∂
x
∂
y
4
+
∂
5
u
∂
x
∂
z
4
)
+
2
∂
5
u
∂
x
∂
y
2
∂
z
2
+
2
l
2
(
∂
3
u
∂
x
∂
y
2
+
∂
3
u
∂
x
∂
z
2
)
+
l
4
(
∂
u
∂
x
)
-
l
2
∂
f
∂
x
+
∂
3
f
∂
x
3
-
∂
3
f
∂
x
∂
y
2
-
∂
3
f
∂
x
∂
z
2
.
Now using center difference discretization at (
0
,
j
,
k
), (25) will become
(65)
g
-
h
2
6
[
∂
2
g
∂
y
2
+
∂
2
g
∂
z
2
+
l
2
(
g
)
-
∂
f
∂
x
]
+
h
4
120
×
[
∂
4
g
∂
y
4
+
∂
4
g
∂
z
4
+
2
∂
4
g
∂
y
2
∂
z
2
+
2
l
2
(
∂
2
g
∂
y
2
+
∂
2
g
∂
z
2
)
+
l
4
(
g
)
-
l
2
∂
f
∂
x
+
∂
3
f
∂
x
3
-
∂
3
f
∂
x
∂
y
2
-
∂
3
f
∂
x
∂
z
2
]
.
Equations (60) and (65) yield that for Neumann boundary points the stencil is of five points for two-dimensional case while in three-dimensional case the stencil has seven points. The matrix is inverted including the extra line i
=
-
1. Applying Neumann boundary conditions accuracy and CPU timing remain the same in all examples.Different methods are used to solve (48) such as finite-difference scheme, fourth-order compact finite-difference scheme, and LU decomposition. We are solving this equation by Multigrid method with sixth-order compact finite-difference scheme using Gauss-Seidel method as a smoother.


## 4. Multigrid Method

The results of second-order finite-difference scheme and sixth-order compact finite-difference scheme are in sparse linear systems; that can be solved efficiently by Multigrid methods. In order to remove high frequency error, Multigrid method utilizes some relaxation methods. Multigrid method makes use of coarse grid correction to smooth the errors. Efficient Multigrid method is implemented in [9, 10] for solving 2D Poisson equation discretized by the standard fourth-order compact difference schemes. It is shown that the sixth-order compact difference scheme is the most cost effective compared to the second-order finite-difference schemes and the corresponding conventional second-order central-difference scheme with Multigrid methods. The standard Multigrid method with a point Gauss-Seidel type relaxation and standard mesh coarsening strategy does not work very well with unequal mesh sizes discretized Poisson equation [11].The strategy used is the line relaxation to replace point relaxation. That is, line Gauss-Seidel relaxation can be shown to be very effective in removing high frequency errors. Consider in this paper for the particular problem the dominant direction is always either thex or y direction, but not both. Let x be the dominant direction; we only perform line relaxation along x direction on each successive grid. The coefficient matrix of the sixth-order compact finite difference with this ordering can be written in block tridiagonal matrix of the block order (
N
-
1
). The order of the coefficient matrix U is of (
N
-
1
)
×
(
N
-
1
), where U
=
diag
⁡
[
U
1
,
U
0
,
U
1
], where U
0
=
diag
⁡
[
a
2
,
a
1
,
a
2
] and U
1
=
diag
⁡
[
a
3
,
a
2
,
a
3
] are both symmetric tridiagonal submatrices of order (
N
-
1
). They represent the submatrix of of each grid line along one direction, where u
j is the part of the solution vector representing the grid points on each jth line, and f
j is the corresponding part on the right-hand side vector. On each level U
0 needs only one factorization. The factorization cost of U
0 is negligible because matrix B has constant blocks which does not change from one grid line to another grid line.In Multigrid method with LU-decompositon or Gauss-Seidel relaxations, we use bilinear interpolation to transfer correction from a coarse grid to a fine grid, and we also use a full-weighting scheme to update the residual on a coarse grid.Algorithm 1 (Multigrid Algorithm).
Assuming that we set up these Multigrid parameters:v
1
pre smoothing steps on each level.
v
2 post smoothing steps on each level.
r is the number of Multigrid cycles on each level (r
=
1 for V-cycling and r
=
2 for W-cycling), here we use V-cycle with r
=
1.
FAS Multigrid Cycle
(66)
ϕ
h
⟵
FASCYC
h
(
ϕ
h
,
f
h
,
v
1
,
v
2
,
r
)
(1)
ifΩ
h is the coarsest grid, then solve (48) using a time marching technique and then stop.
Else do the pre-smoothing step:(67)
ϕ
h
⟵
Smoother
(
ϕ
h
,
f
h
,
v
1
,
tol
)
,
(
Pre-Smoothig
)
.
(2)
Restriction: consider
(68)
ϕ
2
h
=
I
h
2
h
ϕ
h
,
ϕ
2
h
=
ϕ
-
2
h
,
f
2
h
=
I
h
2
h
(
f
h
-
N
h
ϕ
h
)
+
N
2
h
ϕ
2
h
,
ϕ
2
h
⟵
FASCYC
2
h
(
ϕ
2
h
,
f
2
h
,
v
1
,
v
2
)
.
(3)
Interpolation: consider
(69)
ϕ
h
⟵
ϕ
h
+
I
2
h
h
(
ϕ
2
h
-
ϕ
-
2
h
)
.
(4)
One has(70)
ϕ
h
⟵
Smoother
v
2
(
ϕ
h
,
f
h
,
v
2
)
.
(
Post-Smoothing
)
.
Here the restriction operatorI
h
2
h is by full weighting and the interpolation I
2
h
h is by bilinear operator.


## 5. Numerical Calculations

In order to show the efficiency and applicability of the Multigrid method, numerical experiments are conducted to solve a two-dimensional Helmholtz equation (1) on the unit square domain [
0,1
]
×
[
0,1
]. The right hand-side function and the pure Dirichlet boundary conditions are applied on all side of a unit square and unit cubic domain.


### 5.1. 2D Case

Example 1.
Consider the following:(71)
u
x
x
+
u
y
y
+
l
2
u
=
(
l
2
-
2
π
2
)
sin
(
π
x
)
sin
(
π
y
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
.The exact solution isu
(
x
,
y
)
=
sin
(
π
x
)
sin
(
π
y
). The sixth-order compact finite-difference scheme is compared against the standard second-order central difference scheme and fourth-order compact difference scheme, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. One pre-smoothing and post-smoothing are applied at each level. The iteration stops when the Euclidean norm (2-norm) of the residual vector is reduced by 1
0
-
14. The maximum absolute error reported is the maximum absolute error between the computed solution and the exact solution over the entire fine grid points. In order to compare the numerical solution and the exact solution, we use l
2-norm. The matrix l
2-norm of the error vector is defined as follows:
(72)
∥
e
2
∥
=
1
N
∑
i
,
j
=
0
N
‍
e
i
,
j
2
,
where the error vector e
i
j
=
u
i
j
-
v
i
j and the residual r
=
(
f
-
A
v
)
=
A
e, N is the number of nodes, l is the wave number, and furthermore M
1 is the MG with 2nd order, M
2 is MG with 4th order, and M
3 is MG with 6th order. The l
2-norm of the error for l
=
10 and the data in Table 1, Table 2 indicates the behavior of N for different values (see Figure 1).Table 1
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example1, ∥
e
2
∥, l
=
10, and N
=
4,8
,
16,32,64,128.
N
M
1
(
∥
e
2
∥
)
CPU
V
M
2
(
∥
e
2
∥
)
CPU
V
M
3
(
∥
e
2
∥
)
CPU
V
4
1.0000
e
0
0.105
2
2.2284
e
-
4
0.042
2
6.2252
e
-
7
0.042
2
8
3.1345
e
-
3
0.140
2
9.7200
e
-
5
0.055
2
3.2845
e
-
7
0.054
2
16
7.8851
e
-
4
0.780
2
1.9170
e
-
6
0.061
2
4.5025
e
-
9
0.060
2
32
1.9741
e
-
4
0.810
2
6.9331
e
-
8
0.088
2
6.6001
e
-
11
0.087
2
64
4.9377
e
-
5
0.940
2
3.5772
e
-
9
0.322
2
9.9922
e
-
13
0.313
2
128
1.2346
e
-
5
1.137
2
2.1063
e
-
10
2.233
2
1.5369
e
-
14
2.224
2Table 2
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example2, ∥
e
2
∥, l
=
100, and N
=
4,8
,
16,32,64,128.
N
M
1
(
∥
e
2
∥
)
CPU
V
M
2
(
∥
e
2
∥
)
CPU
V
M
3
(
∥
e
2
∥
)
CPU
V
4
6.0562
e
-
2
0.062
2
6.3220
e
-
6
0.042
2
5.6567
e
-
7
0.042
2
8
1.2872
e
-
4
0.074
2
8.1311
e
-
7
0.059
2
7.5902
e
-
8
0.057
2
16
3.2304
e
-
5
0.080
2
3.2134
e
-
7
0.080
2
6.2834
e
-
8
0.063
2
32
8.1531
e
-
6
0.093
2
9.8314
e
-
8
0.187
2
5.5698
e
-
8
0.181
2
64
2.0389
e
-
6
0.118
2
7.0389
e
-
8
0.329
2
5.2102
e
-
8
0.313
2
128
5.0970
e
-
7
0.290
2
6.9061
e
-
8
2.269
2
5.1442
e
-
8
2.252
2(a) The error graph. (b) Exact solution. (c) The approximation solution. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
128 are the number of nodes and l
=
10 for Example 1.
(a)
(b)
(c)Example 2.
Consider the following:(73)
u
x
x
+
u
y
y
+
l
2
u
=
-
2
[
(
1
-
6
x
2
)
(
y
2
-
y
4
)
+
(
1
-
6
y
2
)
(
x
2
-
x
4
)
]
+
l
2
(
(
x
2
-
x
4
)
(
y
4
-
y
2
)
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
.With the pure Dirichlet boundary conditions on all sides of a unit square.The exact solution, for Example2, is u
(
x
,
y
)
=
(
x
2
-
x
4
)
(
y
4
-
y
2
) (see Figure 2).(a) The error graph. (b) Exact solution. (c) The approximation solution. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
128 are the number of nodes and l
=
100 for Example 2.
(a)
(b)
(c)


### 5.2. 3D Case

Example 3.
Consider the following:(74)
u
x
x
+
u
y
y
+
u
z
z
+
l
2
u
(
x
,
y
)
=
(
l
2
-
3
π
2
)
(
cos
⁡
(
π
x
)
sin
(
π
y
)
sin
(
π
z
)
)
,
0
≤
x
,
y
,
z
≤
1
,
with the Neumann boundary condition on the left-hand side of the cubic domain and Dirichlet boundary conditions on the remaining five sides, that is, u
x
(
0
,
y
,
z
)
=
0, u
(
1
,
y
,
z
)
=
-
sin
(
π
y
)
sin
(
π
z
), and u
(
x
,
0
,
z
)
=
u
(
x
,
1
,
z
)
=
u
(
x
,
y
,
0
)
=
u
(
x
,
y
,
1
)
=
0. The exact solution for above example is u
(
x
,
y
,
z
)
=
cos
⁡
(
π
x
)
sin
(
π
y
)
sin
(
π
z
). The Multigrid with sixth-order compact finite-difference scheme is compared against the second-order finite-difference scheme and sixth-order compact finite-difference scheme without Multigrid, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. The iteration stops when the Euclidean norm (2-norm) of the residual vector is reduced by 1
0
-
5. The maximum absolute error is reported in (72). M
4 is the second order finite difference scheme, M
5 is sixth order compact finite difference scheme without Multigrid and M
6 is MG with sixth order compact finite difference scheme. The l
2-norm of the error for l
=
10 and different values of N
=
n-grids, are presented in Table 3.Table 3
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example3, ∥
e
2
∥, l
=
10, and N
=
2,3
,
4,5
,
6.
N
M
4
(
∥
e
2
∥
)
CPU (seconds)
M
5
(
∥
e
2
∥
)
CPU (seconds)
M
6
(
∥
e
2
∥
)
CPU (seconds)
2
8.447
e
-
1
4.130
3.665
e
-
1
2.427
8.454
e
-
3
0.430
3
3.132
e
-
2
6.768
2.762
e
-
2
2.298
3.192
e
-
3
0.677
4
2.152
e
-
2
6.976
9.928
e
-
3
3.010
9.853
e
-
4
0.898
5
8.312
e
-
3
9.203
5.886
e
-
3
3.243
5.229
e
-
4
1.574
6
6.502
e
-
3
14.800
3.990
e
-
3
4.512
1.009
e
-
4
1.906The data in Tables3 and 4 indicates the behavior of N and l for different values in 3-dimensional case. Figures 3 and 4 show the error graph of both problems. Figure 5(a) is the solution of two-dimensional Helmholtz equation while Figure 5(b) is the solution of three-dimensional Helmholtz equation.Table 4
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example4,  ∥
e
2
∥, l
=
10, and N
=
2,3
,
4,5
,
6.
N
M
4
(
∥
e
2
∥
)
CPU (seconds)
M
5
(
∥
e
2
∥
)
CPU (seconds)
M
6
(
∥
e
2
∥
)
CPU (seconds)
2
8.690
e
-
1
3.908
3.000
e
-
3
2.188
1.630
e
-
4
0.039
3
2.429
e
-
1
4.502
5.4781
e
-
4
2.550
8.913
e
-
5
0.317
4
1.109
e
-
2
6.880
2.5575
e
-
4
3.701
5.810
e
-
5
0.856
5
5.616
e
-
3
10.024
1.55553
e
-
4
3.924
2.183
e
-
5
1.538
6
2.700
e
-
3
12.385
1.0792
e
-
4
4.618
9.71
e
-
6
1.876Figure 3
The error graph. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
6 are the number of nodes and l
=
10 for Example 3.Figure 4
The error graph. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
6 are the number of nodes and l
=
10 for Example 4.(a) The approximate solution for 2D. (b) The approximate solution for 3D. ((c) and (d)) The nonzero entries in block tridiagonal matrices, nz are the number of nodes in 2D and 3D cases.
(a)
(b)
(c)
(d)Example 4.
Consider the following:(75)
u
x
x
+
u
y
y
+
u
z
z
+
l
2
u
=
(
1
-
l
2
)
(
sin
(
k
x
2
π
)
sin
(
k
y
2
π
)
sin
(
k
z
2
π
)
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
,
0
≤
z
≤
1
,
with the pure Dirichlet boundary conditions on all sides of a unit square.The exact solution for Example4 is u
(
x
,
y
,
z
)
=
(
sin
(
k
x
2
π
)
sin
(
k
y
2
π
)
sin
(
k
z
2
π
)
)
/
(
k
π
2
(
x
2
+
y
2
+
z
2
)
).


## 5.1. 2D Case

Example 1.
Consider the following:(71)
u
x
x
+
u
y
y
+
l
2
u
=
(
l
2
-
2
π
2
)
sin
(
π
x
)
sin
(
π
y
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
.The exact solution isu
(
x
,
y
)
=
sin
(
π
x
)
sin
(
π
y
). The sixth-order compact finite-difference scheme is compared against the standard second-order central difference scheme and fourth-order compact difference scheme, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. One pre-smoothing and post-smoothing are applied at each level. The iteration stops when the Euclidean norm (2-norm) of the residual vector is reduced by 1
0
-
14. The maximum absolute error reported is the maximum absolute error between the computed solution and the exact solution over the entire fine grid points. In order to compare the numerical solution and the exact solution, we use l
2-norm. The matrix l
2-norm of the error vector is defined as follows:
(72)
∥
e
2
∥
=
1
N
∑
i
,
j
=
0
N
‍
e
i
,
j
2
,
where the error vector e
i
j
=
u
i
j
-
v
i
j and the residual r
=
(
f
-
A
v
)
=
A
e, N is the number of nodes, l is the wave number, and furthermore M
1 is the MG with 2nd order, M
2 is MG with 4th order, and M
3 is MG with 6th order. The l
2-norm of the error for l
=
10 and the data in Table 1, Table 2 indicates the behavior of N for different values (see Figure 1).Table 1
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example1, ∥
e
2
∥, l
=
10, and N
=
4,8
,
16,32,64,128.
N
M
1
(
∥
e
2
∥
)
CPU
V
M
2
(
∥
e
2
∥
)
CPU
V
M
3
(
∥
e
2
∥
)
CPU
V
4
1.0000
e
0
0.105
2
2.2284
e
-
4
0.042
2
6.2252
e
-
7
0.042
2
8
3.1345
e
-
3
0.140
2
9.7200
e
-
5
0.055
2
3.2845
e
-
7
0.054
2
16
7.8851
e
-
4
0.780
2
1.9170
e
-
6
0.061
2
4.5025
e
-
9
0.060
2
32
1.9741
e
-
4
0.810
2
6.9331
e
-
8
0.088
2
6.6001
e
-
11
0.087
2
64
4.9377
e
-
5
0.940
2
3.5772
e
-
9
0.322
2
9.9922
e
-
13
0.313
2
128
1.2346
e
-
5
1.137
2
2.1063
e
-
10
2.233
2
1.5369
e
-
14
2.224
2Table 2
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example2, ∥
e
2
∥, l
=
100, and N
=
4,8
,
16,32,64,128.
N
M
1
(
∥
e
2
∥
)
CPU
V
M
2
(
∥
e
2
∥
)
CPU
V
M
3
(
∥
e
2
∥
)
CPU
V
4
6.0562
e
-
2
0.062
2
6.3220
e
-
6
0.042
2
5.6567
e
-
7
0.042
2
8
1.2872
e
-
4
0.074
2
8.1311
e
-
7
0.059
2
7.5902
e
-
8
0.057
2
16
3.2304
e
-
5
0.080
2
3.2134
e
-
7
0.080
2
6.2834
e
-
8
0.063
2
32
8.1531
e
-
6
0.093
2
9.8314
e
-
8
0.187
2
5.5698
e
-
8
0.181
2
64
2.0389
e
-
6
0.118
2
7.0389
e
-
8
0.329
2
5.2102
e
-
8
0.313
2
128
5.0970
e
-
7
0.290
2
6.9061
e
-
8
2.269
2
5.1442
e
-
8
2.252
2(a) The error graph. (b) Exact solution. (c) The approximation solution. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
128 are the number of nodes and l
=
10 for Example 1.
(a)
(b)
(c)Example 2.
Consider the following:(73)
u
x
x
+
u
y
y
+
l
2
u
=
-
2
[
(
1
-
6
x
2
)
(
y
2
-
y
4
)
+
(
1
-
6
y
2
)
(
x
2
-
x
4
)
]
+
l
2
(
(
x
2
-
x
4
)
(
y
4
-
y
2
)
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
.With the pure Dirichlet boundary conditions on all sides of a unit square.The exact solution, for Example2, is u
(
x
,
y
)
=
(
x
2
-
x
4
)
(
y
4
-
y
2
) (see Figure 2).(a) The error graph. (b) Exact solution. (c) The approximation solution. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
128 are the number of nodes and l
=
100 for Example 2.
(a)
(b)
(c)


## 5.2. 3D Case

Example 3.
Consider the following:(74)
u
x
x
+
u
y
y
+
u
z
z
+
l
2
u
(
x
,
y
)
=
(
l
2
-
3
π
2
)
(
cos
⁡
(
π
x
)
sin
(
π
y
)
sin
(
π
z
)
)
,
0
≤
x
,
y
,
z
≤
1
,
with the Neumann boundary condition on the left-hand side of the cubic domain and Dirichlet boundary conditions on the remaining five sides, that is, u
x
(
0
,
y
,
z
)
=
0, u
(
1
,
y
,
z
)
=
-
sin
(
π
y
)
sin
(
π
z
), and u
(
x
,
0
,
z
)
=
u
(
x
,
1
,
z
)
=
u
(
x
,
y
,
0
)
=
u
(
x
,
y
,
1
)
=
0. The exact solution for above example is u
(
x
,
y
,
z
)
=
cos
⁡
(
π
x
)
sin
(
π
y
)
sin
(
π
z
). The Multigrid with sixth-order compact finite-difference scheme is compared against the second-order finite-difference scheme and sixth-order compact finite-difference scheme without Multigrid, in terms of solution of accuracy, Multigrid convergence rate, and CPU timing. The iteration stops when the Euclidean norm (2-norm) of the residual vector is reduced by 1
0
-
5. The maximum absolute error is reported in (72). M
4 is the second order finite difference scheme, M
5 is sixth order compact finite difference scheme without Multigrid and M
6 is MG with sixth order compact finite difference scheme. The l
2-norm of the error for l
=
10 and different values of N
=
n-grids, are presented in Table 3.Table 3
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example3, ∥
e
2
∥, l
=
10, and N
=
2,3
,
4,5
,
6.
N
M
4
(
∥
e
2
∥
)
CPU (seconds)
M
5
(
∥
e
2
∥
)
CPU (seconds)
M
6
(
∥
e
2
∥
)
CPU (seconds)
2
8.447
e
-
1
4.130
3.665
e
-
1
2.427
8.454
e
-
3
0.430
3
3.132
e
-
2
6.768
2.762
e
-
2
2.298
3.192
e
-
3
0.677
4
2.152
e
-
2
6.976
9.928
e
-
3
3.010
9.853
e
-
4
0.898
5
8.312
e
-
3
9.203
5.886
e
-
3
3.243
5.229
e
-
4
1.574
6
6.502
e
-
3
14.800
3.990
e
-
3
4.512
1.009
e
-
4
1.906The data in Tables3 and 4 indicates the behavior of N and l for different values in 3-dimensional case. Figures 3 and 4 show the error graph of both problems. Figure 5(a) is the solution of two-dimensional Helmholtz equation while Figure 5(b) is the solution of three-dimensional Helmholtz equation.Table 4
Comparison of maximum absolute errors and CPU (seconds) for a Multigrid method with different schemes, for Example4,  ∥
e
2
∥, l
=
10, and N
=
2,3
,
4,5
,
6.
N
M
4
(
∥
e
2
∥
)
CPU (seconds)
M
5
(
∥
e
2
∥
)
CPU (seconds)
M
6
(
∥
e
2
∥
)
CPU (seconds)
2
8.690
e
-
1
3.908
3.000
e
-
3
2.188
1.630
e
-
4
0.039
3
2.429
e
-
1
4.502
5.4781
e
-
4
2.550
8.913
e
-
5
0.317
4
1.109
e
-
2
6.880
2.5575
e
-
4
3.701
5.810
e
-
5
0.856
5
5.616
e
-
3
10.024
1.55553
e
-
4
3.924
2.183
e
-
5
1.538
6
2.700
e
-
3
12.385
1.0792
e
-
4
4.618
9.71
e
-
6
1.876Figure 3
The error graph. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
6 are the number of nodes and l
=
10 for Example 3.Figure 4
The error graph. The error vectore
i
j
=
u
i
j
-
v
i
j and N
=
6 are the number of nodes and l
=
10 for Example 4.(a) The approximate solution for 2D. (b) The approximate solution for 3D. ((c) and (d)) The nonzero entries in block tridiagonal matrices, nz are the number of nodes in 2D and 3D cases.
(a)
(b)
(c)
(d)Example 4.
Consider the following:(75)
u
x
x
+
u
y
y
+
u
z
z
+
l
2
u
=
(
1
-
l
2
)
(
sin
(
k
x
2
π
)
sin
(
k
y
2
π
)
sin
(
k
z
2
π
)
)
,
0
≤
x
≤
1
,
0
≤
y
≤
1
,
0
≤
z
≤
1
,
with the pure Dirichlet boundary conditions on all sides of a unit square.The exact solution for Example4 is u
(
x
,
y
,
z
)
=
(
sin
(
k
x
2
π
)
sin
(
k
y
2
π
)
sin
(
k
z
2
π
)
)
/
(
k
π
2
(
x
2
+
y
2
+
z
2
)
).


## 6. Conclusion

We have studied a sixth-order compact finite-difference scheme with equal mesh sizes for discretizing a two and three-Dimensional Helmholtz equation. We have developed special Multigrid method to solve the resulting sparse system efficiently. It is observed that Multigrid method with the Gauss-Seidel relaxation works very well in solving the sixth-order compact finite-difference scheme-discretized 2D and 3D Helmholtz equation. Numerical results show that Multigrid method on sixth-order compact scheme has the expected accuracy and is faster than the second-order finite-difference scheme and sixth-order compact finite-difference scheme without Multigrid. It is also obvious from the results that the overall error does not decrease with the increase of the wave numberl.



---

*Source: 954658-2014-08-27.xml*