A new 9-point sixth-order accurate

compact finite difference method

for the Helmholtz equation

Majid Nabavi, M. H. Kamran Siddiqui ∗ , Javad Dargahi

Department of Mechanical and Industrial Engineering, Concordia University,
1455 de Maisonneuve Blvd. West, Montreal, Quebec, Canada H3G 1M8

Abstract

A new 9-point sixth-order accurate compact finite difference method for solving the
Helmholtz equation in one and two dimensions, is developed and analyzed. This
scheme is based on sixth-order approximation to the derivative calculated from
the Helmholtz equation. A sixth-order accurate symmetrical representation for the
Neumann boundary condition was also developed. The efficiency and accuracy of
the scheme is validated by its application to two test problems which have exact
solutions. Numerical results show that this sixth-order scheme has the expected
accuracy and behaves robustly with respect to the wave number.

Key words: Helmholtz equation, Compact finite difference, Numerical Analysis

1 Introduction

The Helmholtz equation, (∆ + k 2 )u = f , is an elliptic partial differential equa-

tion which is a time-harmonic solution of the wave equation. The Helmholtz
equation governs some important physical phenomena. These include the po-
tential in time harmonic acoustic and electromagnetic fields, acoustic wave
scattering, noise reduction in silencers, water wave propagation, membrane
vibration and radar scattering. Obtaining an efficient and more accurate nu-
merical solution for the Helmholtz equation has been the subject of many
studies. The numerical solution of the Helmholtz equation has been developed
using different approaches such as the finite difference method [1], the bound-
ary element method [2], the finite element method [3] and the spectral element
∗ Corresponding author
E-mail address: siddiqui@encs.concordia.ca, Tel: 514-8482424 ext. 7940, Fax: 514-8483175

Preprint submitted to Elsevier Science 9 July 2007

method [4].
The boundary element method is derived through the discretization of an
integral equation that is mathematically equivalent to the original partial dif-
ferential equation. The disadvantages of boundary element methods are the
restriction to linear problems in homogeneous and isotropic media, as well as
the large computer storage space required and lengthy processing time needed
to solve the inherent problems encountered with characteristic wave numbers.
Finite element methods are used extensively to solve the Helmholtz equation.
In addition to the high computational cost, another disadvantage of Galerkin
finite element method for solving the Helmholtz equation is the so-called pollu-
tion effect, which results in the less accurate solution at higher wave numbers
for the given nodes per wavelength. Thus, in order to obtain the same solu-
tion accuracy for higher wave numbers, more nodes per wavelength are needed
than that for the lower wave numbers [5, 6]. Although some modification in
the standard Galerkin approximation have been developed to minimize the
pollution effect [8], finding an optimal method is still a challenge [7].
For spectral element methods, it is shown that it requires fewer grid nodes
per wavelength compared to the finite element methods for the Helmholtz
equation [4]. However, due to the less sparse resultant matrix compared to the
resulting finite element matrix, the computational time of both methods are
almost the same [4].
For the traditional finite difference methods, in order to increase the order
of accuracy of approximation, the stencil of grid points needs to be enlarged,
which is not desirable.
Generally, obtaining a more accurate numerical solution means adding more
nodes and using smaller mesh sizes, which requires more computing time and
storage space. In order to obtain more accurate results for constant mesh size,
we need to increase the order of accuracy of the numerical approximation,
which, in turn, means enlarging the stencil of grid points. This, however, leads
to some problems including difficult treatments of the boundary conditions
and approximation of the points next to the boundaries, and increasing the
bandwidth of the stiffness matrix, which makes fast direct solver difficult.
Therefore, compact finite difference schemes are desired to solve partial dif-
ferential equation numerically.
A noticeable work in this field has been done by Turkel and Singer [1]. They
developed a fourth-order compact finite difference method using two schemes.
One scheme was based on the generalization of the Padé approximation and
the other used the Helmholtz equation to calculate higher-order correction
terms. They implemented these schemes for Dirichlet and/or Neumann bound-
ary conditions. In the present study we extended the previous work and devel-
oped a new scheme to increase the accuracy to the order of six without enlarg-
ing the stencil of grid points. We have developed and implemented a new 9-
point stencil, sixth-order accurate compact finite difference method for solving
the Helmholtz equation in the one-dimensional and two-dimensional domain
with Dirichlet and/or Neumann boundary conditions. Recently, Sutmann [9]

2
has reported a sixth-order compact finite difference scheme for Helmholtz
equation. However, he implemented his scheme only with Dirichlet bound-
ary conditions. As far as the authors know, the present study is the first that
developed a sixth-order compact finite difference scheme with the Neumann
boundary conditions.

2 Nine-point sixth-order accurate compact finite difference scheme

This method is based on a sixth-order accurate approximation to the deriva-

tive calculated from the Helmholtz equation. We developed the scheme for
both one-dimensional and two-dimensional uniform Cartesian grids with grid
spacing ∆x = ∆y = h. The notation δx0 is used to denote the first order
central difference with respect to x, which is defined as
ui+1 − ui−1
δx0 ui = , (1)
2h
where ui = u(xi ). The standard second-order central difference is denoted by
δx2 and is defined as
ui+1 − 2ui + ui−1
δx2 ui = . (2)
h2
Difference operators δy0 and δy2 are defined similarly.

2.1 One dimensional case

In one dimensional case, the Helmholtz equation becomes the following ordi-
nary differential equation

u00 + k 2 u = f (x) for x ∈ [a, b]. (3)

A sixth-order accurate finite difference estimate of Eq. (2) is

h2 (4) h4 (6)
δx2 ui = u00i + ui + u + O(h6 ). (4)
12 360 i
Both O(h2 ) and O(h4 ) terms are included in Eq. (4), because we want to
approximate all of them in order to construct an O(h6 ) scheme. Applying δx2
(4)
to ui , we get
(6) (4)
ui = δx2 ui + O(h2 ). (5)
Substituting Eq. (5) into Eq. (4) yields

h2 (4) h4 2 (4)
δx2 ui = u00i + ui + (δ u + O(h2 )) + O(h6 ). (6)
12 360 x i

3
To get a compact O(h6 ) approximation, we simply take the appropriate deriva-
tive of Eq. (3), that is
(4)
ui = −k 2 u00i + fi00 , (7)
where fi = f (xi ) and fi00 = f 00 (xi ). Inserting Eq. (7) into Eq. (6) will result in

³ h2 h4 2 ´
δx2 ui = u00i + + δ (−k 2 u00i + fi00 ) + O(h6 ). (8)
12 360 x

Then, a compact (implicit) approximation for u00i with a sixth-order accuracy

will be given as
³ ´
h2 2
δx2 ui − 12
1 + h30 δx2 fi00
u00i = ³ k2 h2
³
h2 2
´´ + O(h6 ). (9)
1− 12
1 + δ
30 x

Using this estimate and considering the discrete solution of Eq. (3) which
satisfies the approximation, we get

³ k 2 h2 ³ h2 ´´ ³ k 2 h2 ³ h2 ´´ h2 ³ h2 ´
δx2 Ui + k 2 1 − 1 + δx2 Ui = 1 − 1 + δx2 fi + 1 + δx2 fi00 ,
12 30 12 30 12 30
(10)
³ k 2 h2 ´ ³ k 4 h4 ´ 2 ³ k 2 h2 ´ k 2 h4 2 h2 h4 2 00
k 2 1− Ui + 1− δx Ui = 1− fi − δx fi + fi00 + δ f ,
12 360 12 360 12 360 x i
(11)
where Ui is the discrete approximation to ui satisfying the discrete formulation
of Eq. (3) which implies, ui = Ui + O(h6 ). Using Eq. (2) and δx2 fi00 = (fi+1 00
−
00 00 2
2fi + fi−1 )/h , we can express the scheme in the following form

d10 Ui +d11 (Ui+1 +Ui−1 ) = b10 fi +b11 (fi+1 +fi−1 )+b12 fi00 +b13 (fi+1
00 00
+fi−1 ), (12)

where,

³28k 2 h2 ´ k 4 h4
d10 = −2 + k 2 h2 1 − , d11 = 1 − ,
360 360
³ 28k 2 h2 ´ 2 −k 2 h4 28h4 h4
b10 = 1− h , b11 = , b12 = , b13 = . (13)
360 360 360 360

Since f and f 00 are known at every grid point, the right hand side of Eq.
(12) is known for all nodes. The system equations given by Eq. (12) can be
written for each node and a resultant linear system of equation is obtained.
In the cases where f is not known analytically, only a fourth-order accurate
approximation of f 00 is required, which can be obtaind using fi00 = (−fi+1 +
16fi+1/2 − 30fi + 16fi−1/2 − fi−1 )/12h2 .

4
2.2 Two dimensional case

Consider the two-dimensional Helmholtz equation

∂ 2u ∂ 2u
+ + k 2 u(x, y) = f (x, y) for x ∈ [a, b] and y ∈ [c, d]. (14)
∂x2 ∂y 2
We want to have a 9-point stencil which is symmetrical in both x and y
directions. The central difference scheme for Eq. (14) in two dimensions can
be written as
δx2 ui,j + δy2 ui,j + k 2 ui,j + Ti,j = fi,j , (15)
where ui,j = u(xi , yj ), fi,j = f (xi , yj ) and

h2 h ∂ 4 u ∂ 4 u i h4 h ∂ 6 u ∂ 6 u i
Ti,j = − + − + + O(h6 ). (16)
12 ∂x4 ∂y 4 i,j 360 ∂x6 ∂y 6 i,j
Using the following appropriate derivatives of Eq. (14)

∂ 4u ∂2f 2
2∂ u ∂4u ∂ 4u ∂ 2f 2
2∂ u ∂ 4u
= − k − , = − k − (17)
∂x4 ∂x2 ∂x2 ∂x2 ∂y 2 ∂y 4 ∂y 2 ∂y 2 ∂y 2 ∂x2
in Eq. (16), we get

h2 ³ 2 h ∂ 4u i h 2
2 ∂ u ∂ u
2 i ´
h4 h ∂ 6 u ∂ 6 u i
Ti,j = − ∇ fi,j −2 2 2
−k 2
+ 2 − 6
+ 6 +O(h6 ).
12 ∂x ∂y i,j ∂x ∂y i,j 360 ∂x ∂y i,j
(18)
In our derivation of an O(h6 ) scheme, we need a fourth-order approximation
of ∂ 4 u/∂x2 ∂y 2 in Eq. (18) which can be written as
h ∂ 4u i 2 2 h2 h ∂ 6 u ∂ 6u i
= δ δ
x y u i,j − + + O(h4 ). (19)
∂x2 ∂y 2 i,j 12 ∂x4 ∂y 2 ∂x2 ∂y 4 i,j
Substituting Eq. (19) into Eq. (18), we get

h2 ³ ´
Ti,j = − ∇2 fi,j + 2δx2 δy2 ui,j + k 2 fi,j − k 4 ui,j
12
h4 h ∂ 6 u ∂6u ∂ 6u ∂ 6u i
− 6
+ 5 4 2
+ 5 2 4
+ 6
+ O(h6 ). (20)
360 ∂x ∂x ∂y ∂x ∂y ∂y i,j

Clearly, getting a compact sixth-order approximation requires compact ex-

pressions of the four derivatives of order six in Eq. (20), which can be done
by further differentiating Eq. (14), that is

∂ 4f ∂ 6u ∂ 6u 4
2 ∂ u
= + + k , (21)
∂x2 ∂y 2 ∂x4 ∂y 2 ∂x2 ∂y 4 ∂x2 ∂y 2
∂ 6u ∂ 6u ³ 4 ∂ 4u ´ ³ ∂ 6u ∂ 6u ´
4 2 ∂ u
+ = ∇ f − k + − + . (22)
∂x6 ∂y 6 ∂x4 ∂y 4 ∂x4 ∂y 2 ∂x2 ∂y 4

5
Substituting Eqs. (17,21) into Eq. (22) gives us

∂ 6u ∂6u 4 ∂ 4f 2 2 4 2
4
2 ∂ u
+ = ∇ f − − k ∇ f + k (−k u + f ) + 3k . (23)
∂x6 ∂y 6 ∂x2 ∂y 2 ∂x2 ∂y 2
Using Eqs. (21,23), we can eliminate all the derivatives of u in Eq. (20), that
is
h2 ³ ´
Ti,j = − ∇2 fi,j + 2δx2 δy2 ui,j + k 2 fi,j − k 4 ui,j
12
h4 ³ 4 h ∂ 4f i
2 2 4 6 2 2 2
´
− ∇ fi,j +4 −k ∇ fi,j +k f i,j −k u i,j −2k δ δ
x y u i,j +O(h6 ).
360 ∂x2 ∂y 2 i,j
(24)

The compact sixth-order approximation of the two-dimensional Helmholtz

equation can thus be obtained as:

h2 ³ k 2 h2 ´ 2 2 ³ k 2 h2 k 4 h4 ´
1+ δx δy Ui,j + (δx2 + δy2 )Ui,j + k 2 1 − + Ui,j =
6 30 12 360
³ k 2 h2 k 4 h4 ´ ³ h2 ³ k 2 h2 ´´ 2 h4 4 h4 h ∂ 4 f i
1− + fi,j + 1− ∇ fi,j + ∇ fi,j + ,
12 360 12 30 360 90 ∂x2 ∂y 2 i,j
(25)

where Ui,j is the discrete approximation to ui,j satisfying the discrete formula-
tion of Eq. (14) which means ui,j = Ui,j + O(h6 ). As it is seen, we can express
the equation in the form of
h ∂ 4f i
d20 Ui,j + d21 D1 + d22 D2 = b20 fi,j + b21 ∇2 fi,j + b22 ∇4 fi,j + b23 , (26)
∂x2 ∂y 2 i,j

where,

D1 = Ui+1,j + Ui−1,j + Ui,j+1 + Ui,j−1 ,

D2 = Ui+1,j+1 + Ui−1,j+1 + Ui+1,j−1 + Ui−1,j−1 ,
B1 = fi+1,j + fi−1,j + fi,j+1 + fi,j−1 ,
B2 = fi+1,j+1 + fi−1,j+1 + fi+1,j−1 + fi−1,j−1 , (27)

then we get
10 ³ 46 k 2 h2 k 4 h4 ´ 2 k 2 h2 1 k 2 h2
d20 =− + k 2 h2 − + , d21 = − , d22 = + ,
3 45 12 360 3 90 6 180
³ k 2 h2 k 4 h4 ´ h4 ³ k 2 h2 ´ h6 h6
b20 =h2 1 − + , b21 = 1− , b22 = , b23 = . (28)
12 360 12 30 360 90
If we write the system equations given by Eq. (26) for each node, we can
obtain the final linear system of equations. In the cases where f is not known
analytically, only a fourth-order accurate approximation of ∇2 f and a second-
4f
order accurate approximation of ∇4 f and ∂x∂2 ∂y 2 are required.

6
2.2.1 Sixth-order accurate approximation of the boundary

We like to approximate the boundary conditions with the same accuracy as

the interior nodes. For the Dirichlet boundary condition, we can simply put
the boundary for every node on the boundary. For the Neumann boundary
condition, the sixth-order approximation is conducted for both one dimen-
sional and two dimensional cases.

One dimensional case

For a Neumann boundary condition in one dimension, we assume

u0 (x0 ) = β (β a constant). (29)

The sixth-order approximation of Eq. (29) is

h2 000 h4 (5)
δx0 ui = u0i + ui + ui + O(h6 ), (30)
6 120
h2 000 h4 ³ 2 000 ´
δx0 ui = u0i + ui + δx ui + O(h2 ) + O(h6 ). (31)
6 120
Differentiating Eq. (3), we get

u000 2 0 0
i = −k ui + f . (32)

Using Eq. (32) in Eq. (31), we get

k 2 h2 0 k 2 h4 2 0 h2 h2
δx0 ui = (1 − )u − δx ui + (1 + δx2 )fi0 + O(h6 ). (33)
6 120 6 20
Using δx2 u0i = u000 2
i + O(h ) and Eq. (32) in Eq. (33) gives us

k 2 h2 k 4 h4 0 h2 ³ h2 ´ 0 h4 2 0
δx0 ui = (1 − + )ui + 1− fi + δ f + O(h6 ). (34)
6 120 6 20 120 x i
Using δx2 fi0 = (fi+1
0
− 2fi0 + fi−1
0
)/h2 + O(h2 ) and Eq. (1) will result in

ui+1 − ui−1 k 2 h2 k 4 h4 0 h2 ³ 9 k 2 h2 ´ 0 h2 0
= (1− + )u + − fi + (f +f 0 )+O(h6 ).
2h 6 120 6 10 20 120 i+1 i−1
(35)
Considering discrete formulation and using Eq. (29) for i = 0, we get
³k 2 h2 k 4 h4 ´ h3 ³ 9 k 2 h2 ´ 0 h3
d11 (U1 −U−1 ) = 2hβd11 1− + + d11 − f0 + d11 (f10 +f−1
0
).
6 120 3 10 20 60
(36)
We wish to eliminate U−1 because we do not have equations at point x−1 .
Using Eq. (12) for i = 0, we get:

d11 (U1 + U−1 ) + d10 U0 = b10 f0 + b11 (f1 + f−1 ) + b12 f000 + b13 (f100 + f−1
00
) (37)

7
We use Eqs. (36,37) to eliminate U−1 and get the desired approximation for
the boundary point x0 , that is

³
k 2 h2 k 4 h4 ´
2d11 U1 + d10 U0 = 2hβd11 1 − + + b10 f0 + b11 (f1 + f−1 )
6 120
h3 ³ 9 k 2 h2 ´ 0 h3
+ d11 − f0 + d11 (f10 + f−1
0
) + b12 f000 + b13 (f100 + f−1
00
). (38)
3 10 20 60

All parameters on the right hand side of Eq. (38) are known. In the cases where
f is not known analytically, only a fourth-order accurate approximation of f 0
is required.

Two dimensional case

For a Neumann boundary condition in two dimensions, we assume

∂u ¯¯
¯ = β(y). (39)
∂x x=0

The sixth-order approximation of Eq. (39) is

h ∂u i h2 h ∂ 3 u i h4 h ∂ 5 u i
δx0 ui = + + + O(h6 ). (40)
∂x i,j 6 ∂x3 i,j 120 ∂x5 i,j

Using the following appropriate derivatives of Eq. (14), we get

∂ 3u ∂f 2 ∂u ∂ 3u ∂ 5u ∂ 3f 3
2∂ u ∂ 5u
= −k − , = −k − , (41)
∂x3 ∂x ∂x ∂x∂y 2 ∂x5 ∂x3 ∂x3 ∂x3 ∂y 2

In our derivation of an O(h6 ) scheme, we need a fourth-order approximation

of (∂ 3 u/∂x∂y 2 ) in Eq. (41) which can be written as

h ∂ 3u i h2 h ∂ 5 u ∂5u i
= δx δy2 ui,j − + 2 + O(h4 ). (42)
∂x∂y 2 i,j 12 ∂x∂y 4 ∂x3 ∂y 2 i,j

Substituting of Eq. (42) into Eq. (41) gives us

h ∂ 3u i h ∂f i h ∂u i h2 h ∂ 5 u ∂ 5u i
3
= − k2 − δx δy2 ui,j +
4
+ 2 3 2
+ O(h4 ).
∂x i,j ∂x i,j ∂x i,j 12 ∂x∂y ∂x ∂y i,j
(43)
The second-order approximation of (∂ 3 u/∂x3 ) can be written as

h ∂ 3u i h ∂f i h ∂u i
= − k2 − δx δy2 ui,j + O(h2 ). (44)
∂x3 i,j ∂x i,j ∂x i,j

8
Using Eqs. (41,44), the second-order approximation of (∂ 5 u/∂x5 ) can be writ-
ten as
h ∂ 5u i h ∂ 3f i h ∂f i h ∂u i h ∂5u i
= − k2 + k4 + k 2 δx δy2 ui,j − + O(h2 ).
∂x5 i,j ∂x3 i,j ∂x i,j ∂x i,j 3
∂x ∂y 2 i,j
(45)
The appropriate derivatives of Eq. (14) gives us:

∂ 3f ∂ 5u ∂ 5u 3
2 ∂ u ∂ 5u ∂ 5u ∂3f 3
2 ∂ u
= + +k ⇒ + = −k .
∂x∂y 2 ∂x3 ∂y 2 ∂x∂y 4 ∂x∂y 2 ∂x3 ∂y 2 ∂x∂y 4 ∂x∂y 2 ∂x∂y 2
(46)
Substituting Eqs. (43,45,46) into Eq. (40), and after some modifications we
get

k 4 h4 h ∂u i h2 ³ k 2 h2 ´ 2
δx0 ui − + 1+ δx δy ui,j =
120 ∂x i,j 6 30
³ k 2 h2 ´h ∂u i h2 ³ k 2 h2 ´h ∂f i h4 h ∂ 3 f i h4 h ∂ 3 f i
1− + 1− + + +O(h6 ).
6 ∂x i,j 6 20 ∂x i,j 120 ∂x3 i,j 72 ∂x∂y 2 i,j
(47)

Now for (∂u/∂x)i,j on the right hand side of Eq. (47), we use this approxima-
tion
[∂u/∂x]i,j = δx0 ui + µh2 δx δy2 ui,j , (48)
where µ is an arbitrary number. Using Eq. (48) in Eq. (47), multiplying by
2h and setting i = 0 will result in

dˆ21 (U1,j − U−1,j ) + dˆ22 (U1,j+1 + U1,j−1 − U−1,j+1 − U−1,j−1 ) =

³ k 2 h2 ´ h2 ³ k 2 h2 ´h ∂f i h4 h ∂ 3 f i h4 h ∂ 3 f i
1− βj + 1− + + . (49)
6 6 20 ∂x 0,j 120 ∂x3 0,j 72 ∂x∂y 2 0,j

where βj = β(yj ) and

ˆ k 4 h4 1 ³ k 2 h2 ´ ˆ k 4 h4 1 ³ k 2 h2 ´
d21 = 1− (1−2µ)− 1+ , d22 = − µ+ 1+ . (50)
120 6h 30 120 12h 30
Setting i = 0 in Eq. (26), we get

d20 U0,j +d21 (U1,j +U−1,j +U0,j+1 +U0,j−1 )+d22 (U1,j+1 +U−1,j+1 +U1,j−1 +U−1,j−1 )
h ∂ 4f i
= b20 f0,j + b21 ∇2 f0,j + b22 ∇4 f0,j + b23 . (51)
∂x2 ∂y 2 0,j

In order to eliminate all elements with i = −1, we define a constant η such

that η = d21 /dˆ21 = d22 /dˆ22 , which can be obtained by

1
η= 4 h4 . (52)
1 − k120

9
If we multiply Eq. (49) with η and add to Eq. (51) we will get this final formula
for boundary nodes

d20 U0,j + d21 (2U1,j + U0,j+1 + U0,j−1 ) + d22 (U1,j+1 + U1,j−1 ) =

h ∂ 4f i
b20 f0,j + b21 ∇2 f0,j + b22 ∇4 f0,j + b23 +
∂x2 ∂y 2 0,j
³³ k 2 h2 ´ h2 ³ k 2 h2 ´h ∂f i h4 h ∂ 3 f i h4 h ∂ 3 f i ´
η 1− βj + 1− + + .
6 6 20 ∂x 0,j 120 ∂x3 0,j 72 ∂x∂y 2 0,j
(53)

All parameters on the right hand side of Eq. (53) are known. In the cases where
f is not known analytically, only a fourth-order accurate approximation of ∂f∂x
∂3f ∂3f
and a second-order approximation of ∂x3 and ∂x∂y2 are required.

3 Numerical Results

In order to validate our sixth-order accurate scheme and examine its behavior,
we developed the scheme on two model problems in two dimensions. In problem
A, we solved

∂ 2u ∂ 2u 2
2
+ 2 +k u(x, y) = (k 2 −2π 2 ) sin(πx) sin(πy) 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1,
∂x ∂y
(54)
with the pure Dirichlet boundary conditions on all sides of a unit square, that
is
u(0, y) = u(1, y) = u(x, 0) = u(x, 1) = 0, (55)
and in problem B, we solved

∂ 2u ∂2u 2
2
+ 2 +k u(x, y) = (k 2 −2π 2 ) cos(πx) sin(πy) 0 ≤ x ≤ 1 and 0 ≤ y ≤ 1,
∂x ∂y
(56)
with the Neumann boundary condition on the left side of a unit square and
Dirichlet boundary conditions on the remaining three sides, that is

ux (0, y) = 0 u(1, y) = − sin(πy) u(x, 0) = u(x, 1) = 0. (57)

The exact solutions for problems A and B are u(x, y) = sin(πx) sin(πy) and
u(x, y) = cos(πx) sin(πy), respectively.
The next step is to solve the resultant linear set of equations. We used LU-
decomposition by Gaussian elimination with pivoting for each set of equations
a to find values of u at N × N nodes. The code was written in MATLAB
environment using version 7.1. The code was executed on a Pentium 4, 3Ghz
PC.

10
In order to compare the numerical solution Ui,j to the exact solution ui,j , we
use two performance metrics namely l2 -norm and order. The metric l2 -norm
of the error vector e, is defined by
v
u
1u N
uX 2
kek2 = t ei,j , (58)
N i,j=0

where ei,j = ui,j − Ui,j and N is the number of nodes. The metric order is
defined as
kek∞ (n)
Order(n, n + 1) = log2 , (59)
kek∞ (n + 1)
and measures the order of accuracy of numerical solutions. kek∞ in Eq. (59)
is called l∞ -norm of the error vector and is defined as

kek∞ = max ei,j (60)

0≤i,j≤N

The l2 -norm of the error and the order of accuracy, for different values of N and
for k=10, are presented in tables 1 and 2 for problems A and B, respectively. It
is clearly seen that the norm of the error behaves like the order of the scheme.
As we multiply N by two, the error decreases by 26 = 64. It means that our
scheme has really the accuracy of order six. This trend can also be seen in
figs. 1 and 2, where the log2 kek∞ is plotted versus log2 N for both problems
A and B, respectively. The slope of the line is -6 which means that the order
of accuracy is 6.
We also examined the behavior of our scheme for different value of k. Figs.
3 and 4 show log2 kek2 versus k with three different value of N for both
problems A and B, respectively. Fig. 3 shows that for problem A, except for
4 ≤ k ≤ 5 in which the scheme is more sensitive to the value of k, the
method behaves robustly with respect to the wave number. Fig. 4 shows that
compared to problem A, the scheme for problem B is more sensitive to the
value of k. However, for any given value of N , the overall error does not
increase with the value of k. As the figures show, at some particular values
of k in both problems, the accuracy of the method is poor. This fact can
be explained through eigenvalue analysis. The approximate eigenvalues for
2 2 2 2
problem ³ A are given as,´ λAm,n = k − π (m + n ) and for problem B, λBm,n =
k 2 −π 2 (m+1/2)2 +n2 (see [1]). For example, near k = 4.44 where λA1,1 → 0,
problem A is unstable and near k = 5.663 where λB1,1 → 0, problem B is
unstable and the scheme has poor accuracy.
One of the important advantages of the proposed scheme is that in comparison
with the finite element, boundary element or spectral element methods, this
method is very fast. A quantitative comparison in terms of the execution time
between the present scheme and fourth-order accurate finite element method
(FEM) is presented in the last two columns of tables 1 and 2. The results
show that for large number of nodes, the present scheme is more that 100
times faster than the fourth-order FEM.

11
4 Conclusions

A new 9-point sixth-order accurate compact finite difference scheme for the
Helmholtz equation was developed. A sixth-order accurate symmetrical repre-
sentation for the Neumann boundary condition was also developed. Numerical
results show that our scheme has the expected accuracy and is more than 100
times faster than the fourth-order FEM. The results also show that the overall
error does not increase with an increase in the wave number.

Acknowledgements

This research is funded by the grants from Natural Science and Engineering
Research Council of Canada (NSERC) and Concordia University.

References

[1] I. Singer and E. Turkel, High-order finite difference method for the helmholtz
equation, Computer methods in applied mechanics and engineering 163 (1998)
343–358.
[2] R.P. Shaw, Integral equation methods in acoustics, Boundary elements 4 (1988)
221–244.
[3] I. Harari and T.J.R. Hughes, Finite element methods for the helmholtz equation
in an exterior domain: model problem, Computer methods in applied mechanics
and engineering 87 (1991) 59–96.
[4] O. Mehdizadeh and M. Paraschivoiu, Investigation of a two-dimensional spectral
method for helmholtz’s equation, Journal of computational physics 189 (2003)
111–129.
[5] F. Ihlenburg, I. Babuska, Finite element solution to the Helmholtz equation with
high wavenumber part I: the h-version of the FEM, Computers and Mathematics
with applications 30 (1995) 9–37.
[6] F. Ihlenburg, I. Babuska, Finite element solution of the Helmholtz equation
with high wave number part II: the h-p-version of the FEM, Computers and
Mathematics with applications 34 (1997) 315–358.
[7] I. Harari and C.L. Nogueria, Reducing dispersion of linear triangular elements
for the helmholtz equation, Journal of Engineering Mechanics 128 (2002) 351–
358.
[8] A. Oberai and P. Pinsky, A numerical comparison of finite element methods for
the helmholtz equation, Journal of computational Acoustics 8 (2000) 211–221.
[9] G. Sutmann, Compact inite difference schemes of sixth order for the helmholtz
equation, Journal of Computational and Applied Mathematics, (in press).

12
 List of tables

Table 1
Numerical results for the problem A, k = 10
N kek2 kek∞ Order Execution time of the Execution time of the
proposed scheme (sec) 4-order FEM (sec)
8 1.8561E-7 3.7586E-6 0.021 0.84
16 1.4556E-9 5.2585E-8 6.16 0.031 3.66
32 1.1750E-11 7.9975E-10 6.04 0.094 15.86
64 9.4278E-14 1.2433E-11 6.00 0.50 84.14
128 9.5151E-016 1.9691E-013 5.98 4.18 502.84

Table 2
Numerical results for the problem B, k = 10
N kek2 kek∞ Order Execution time of the Execution time of the
proposed scheme (sec) 4-order FEM (sec)
8 2.611E-7 3.368E-6 0.021 0.86
16 2.048E-9 8.112E-8 6.05 0.031 3.71
32 1.653E-11 1.2337E-9 6.04 0.11 17.47
64 1.3264E-13 1.9248E-11 6.00 0.58 95.11
128 1.0560E-015 3.0507E-013 5.98 4.42 541.15

13
 −15

−20

−25
log2(Li−norm)

−30

−35

−40

−45
3 3.5 4 4.5 5 5.5 6 6.5 7
log2(N)

Fig. 1. log2 kek∞ versus log2 N for problem A

−15

−20

−25
log2(Li−norm)

−30

−35

−40

−45
3 3.5 4 4.5 5 5.5 6 6.5 7
log2(N)

Fig. 2. log2 kek∞ versus log2 N for problem B

14
 −20

−25
Log2(L2−norm)

−30

−35

−40

2 4 6 8 10 12 14 16 18 20
k

Fig. 3. log2 kek2 versus k for problem A ; • , N = 8 ; × , N = 16 ; ◦ , N = 32 ; k

varies in units of 0.2.

−15

−20

−25
Log2(L2−norm)

−30

−35

−40

−45
2 4 6 8 10 12 14 16 18 20
k

Fig. 4. log2 kek2 versus k for problem B ; • , N = 8 ; × , N = 16 ; ◦ , N = 32 ; k

varies in units of 0.2.

15