Acoustic Green’s Functions using the 2D

Sinc-Galerkin Method
Adrian R. G. Harwood* and Iain D. J. Dupère

Abstract. In many acoustic problems, the radiated sound field is dominated by scattering effects.
Green’s functions represent the scattering behaviour of a particular geometry and are required to
propagate acoustic disturbances through complex geometries using integral methods. The versatil-
ity of existing integral methods of acoustic propagation may be greatly increased by using numerical
Green’s functions computed for more general geometries. To this end, we investigate the use of the
Sinc-Galerkin method to compute the Green’s function for the Helmholtz equation subject to ho-
mogeneous Dirichlet boundary conditions. We compare the results to a typical boundary element
method implementation. The Sinc-Galerkin procedure demonstrates improved performance on a
number of configurations tested in comparison to the BEM. In particular, accuracy comparable
to BEM can be achieved in far less time while being less sensitive to both frequency and source
position. Although the BEM captures the tip of the singularity more completely, the Sinc-Galerkin
is seen to remain robust despite close proximity of the source point to the boundaries and a coarse,
non-uniform distribution of mesh nodes. However, the characteristic exponential convergence is
slower than many other Sinc-Galerkin applications (or lost entirely) due to the presence of the
domain singularity typical of Green’s functions. In addition, we present one possible formulation
for homogeneous Neumann boundary conditions (rigid boundaries).
Keywords. Sinc Methods, Sinc-Galerkin, Boundary Value Problem, Green’s Functions, Wave prob-
lems.

1. Introduction
The noise associated with engineering products can be a source of annoyance, distraction or indeed a
hazard to health [12]. In particular, motivated by the increase in jet air travel in the 1950s, there has
been a focussed effort on developing mathematical means of describing the mechanism of noise gener-
ation associated with unsteady fluid flow (aero-acoustic noise) [17, 18]. This description together with
a means of describing the propagation and scattering of the noise by solid or fluid boundaries allows
us to design quieter products in the future. One prominent achievement has been the development
of a range of noise prediction schemes which make use of numerical methods, analytical models, or a
combination of both, to calculate both the hydrodynamic and acoustic pressure fields from a source
region all the way to an observer location [47]. A particular class of acoustic prediction schemes,
known as hybrid schemes, use a combination of numerical techniques and analytical formulae to com-
pute an aero-acoustic source and propagate it to an observer, usually many wavelengths away [39].
Use of these schemes as an alternative to solely numerical procedures precludes the need to mesh and
solve equations in the extensive region between the source and the observer meaning the demand on
computational resources is significantly less. However, hybrid schemes generally rely on the evaluation
of an integral equation for the propagation step whose kernel consists of a Green’s function. This

* Corresponding Author.
2 Harwood and Dupère

Green’s function often satisfies the boundary conditions in the propagation domain. Physically, this
Green’s function encapsulates all the complex interactions of the source as the radiation propagates
from a particular source location to a particular observer location. Therefore, the Green’s function is
a powerful representation which, when convoluted with a description of the source, may be used to
yield the observed acoustic field in a single integration. Analytical representations of Green’s functions
are few and far between, and are generally approximations [13]. Hence we are limited to propagating
sound using integral equations in regions where Green’s functions are available. In order to extend the
versatility of existing hybrid CAA schemes, we develop here a procedure suitable for computing the
Green’s function for time harmonic acoustic problems. The inclusion of this procedure in a patching
algorithm will allow the accurate and robust computation acoustic Green’s functions in geometries
where analytical representations are unavailable. Alternatively, the method may be used to compute
Green’s functions for periodic domains. The robustness of current techniques are often hampered by
conditional or slow convergence [14, 19, 20, 29].
1.1. Sinc Methods
The set of Sinc numerical methods represent the function of interest as an expansion of translated
Sinc functions. For the Paley-Wiener class of functions (see Lund and Bowers [21] for their properties),
this infinite series representation is exact and called the Whittaker-Cardinal function. The represen-
tation was initially presented by Whittaker [49] but more recently has been developed into a practical
spatial/spectral method for the solution of elliptic and parabolic linear and non-linear differential
equations by Stenger [43]. This is facilitated by using the Sinc expansion to approximate functions
not of the Paley-Wiener class. This approximate representation then forms the basis of associated nu-
merical methods. Stenger demonstrates that for appropriately selected discretisation parameters, the
unique properties of the Sinc basis ensure exponential convergence of the Sinc series representation.
In one dimension, the Sinc basis functions exist along the whole real line, centred at uniformly dis-
tributed sinc nodes. This infinite interval, and the functions within the space, may be transformed to
any semi-infinite or finite interval through the specification of appropriate 1D conformal maps. This
allows functions on an arbitrary 1D interval, as well as their integrals and derivatives, to be approxi-
mated using a transformed expansion of Sinc functions centred at the transformed nodal locations in
the interval. For many practical problems, such as the ones presented in this paper, the solution in
the domain is assumed unknown. The Sinc expansion may be then used in collocation [32, 33, 50, 51],
Galerkin [7, 26] or highest derivative [40] procedures to find the unknown coefficients. The collocation
method is generally chosen over the Galerkin method for problems with variable coefficients as the
resulting implementation is a more efficient means of achieving a given accuracy [50, 51]. Despite this
benefit, the Galerkin procedure is more flexible as it lends itself to basis modification and extension
without the stringent smoothness requirements of the collocation method [28] and so is used here. A
rigorous mathematical summary of the Sinc method formulae can be found in Lund and Bowers [21]
and Stenger [44].
An important practical advantage of using spectral-like methods over other numerical methods such as
the boundary element method (BEM) is that once the coefficients have been found, the evaluation of
points within and on the domain boundaries requires simply the evaluation and assembly of the series
components. We compare our Sinc implementation to a typical BEM implementation in Section 6,
where the issue of re-evaluating coefficient matrices for every internal point is seen to significantly
impact execution time.
The principal advantage of the Sinc-Galerkin method over Galerkin methods with different basis func-
tions is the ability to achieve exponential convergence of the series solution. This behaviour has been
established for a number of problems in the literature, even for those featuring a boundary singularity
[6, 23, 26, 41]. However, the method has not been applied to a BVP featuring a domain singularity. We
present such an application here. Ensuring exponential convergence means that, compared with other
numerical approaches, a fewer number of functions and nodal points need be used to accurately rep-
resent the function of interest [33]. The single exponential convergence of early Sinc implementations
may by improved upon by using a double exponential conformal mapping [34, 45, 50, 51]. Additional
 Green’s Functions using Sinc-Galerkin 3

advantages of the Sinc-Galerkin approach when compared with alternative numerical methods include
[3]:
1. Problems in semi-infinite and infinite domains may be solved without the need for artificial
boundaries and non-reflecting boundary conditions which are almost always inexact and can be
complicated to implement;
2. Sinc methods do not require numerical integration in order to compute the component matrices,
just point evaluations, which require significantly less computational effort;
3. The error of the approximation is not dependent on any derivatives of the exact solution hence
it is capable of approximating singular solutions without amplification of errors;
4. The exponential convergence rate is maintained even for problems with boundary singularities
[26];
5. Time-dependent problems may be expressed as a linear system rather than a system of ODEs
typical of many Galerkin procedures [26];
6. The weight function in the Galerkin procedure can be selected to ensure a symmetric matrix
system which can increase solution efficiency [25] as well as to [5]:
(a) Simplify the discrete system;
(b) Introduce special structure to reduce computational cost;
(c) Reduce the truncation error of the infinite series;
(d) Force boundary terms to zero (necessary assumption during the application of quadrature);
(e) Allow the approximation to match physical solutions more closely thus reducing overall
errors due to coefficient inaccuracies;
(f) Allow the direct analysis of problems on semi-infinite and infinite domains [32];
However, unlike the Finite Element Method (FEM) and the Finite Difference Method (FDM), the en-
suing matrices are full which precludes the use of efficient sparse solvers. According to McArthur et al.
[25, 26], this disadvantage only partially offsets
√
the significant benefit of the exponential convergence
rate of the Sinc series equivalent to O(e−k N ) where k > 1 and the number of basis functions used to
construct the approximate solution is 2N + 1. In addition, the matrix systems, although not banded,
consist of Toeplitz matrices and, for linear problems, the form of the resulting matrix equation is that
of the Sylvester equation [42]. This particular equation may be solved using a number of established
techniques; two of the most common being described by Lund and Bowers [21], with efficient solvers
available, for handling the special structure of the constituent matrices.
The merits of the Sinc-Galerkin method have led to its applications to both 1D and higher-dimensional
problems subject to homogeneous Dirichlet [3, 5–11, 15, 16, 22, 23, 25, 26, 28, 31–34, 36, 37, 41, 42,
52, 53] and inhomogeneous Dirichlet boundary conditions [7, 15, 25, 31, 34, 41, 42] for a variety of dif-
ferential operators. The wave equation and the Helmholtz equation, previously treated by McArthur,
Bowers and Lund [26] (wave equation) and [3] (Helmholtz as a general elliptic equation), are of par-
ticular interest in acoustics and will be the focus of this paper. The use of Sinc methods with Dirichlet
conditions is relatively straight-forward since the Sinc function satisfies the homogeneous Dirichlet
condition automatically. A simple change of variable is all that is required for the conversion of the
problem with inhomogeneous conditions to the problem with homogeneous conditions as discussed
later. However, when the boundary conditions are Neumann conditions (or mixed conditions involv-
ing both the solution and its derivative), the Sinc basis may be modified to ensure the derivative of the
individual basis functions exist at the endpoints, since ordinarily the derivative of the Sinc function
is undefined here. These modified functions, however, although defined at the boundary, still do not
satisfy the boundary conditions. Therefore, the basis may be extended by including additional basis
functions to interpolate the boundary conditions [21]. Such functions are typically Hermite polynomi-
als. Many examples of the modification [15, 33, 42] and the extension [5, 15, 22, 23, 28, 33, 42] of the
basis to treat a variety of BVPs exist. One apparent absentee is that of the Sinc-Galerkin formulation
for the Helmholtz equation on a rectangle subject to Neumann or mixed conditions on the boundary.
The physical significance of these problems in acoustics are those subject to rigid wall or impedance
boundary conditions. We present one possible formulation in Appendix A. Since this problem is a
4 Harwood and Dupère

‘pure Neumann problem’, it is expected that the compatibility condition must also be satisfied to
ensure that the solution to the problem remains unique; the boundary conditions not being sufficient
to fully define the solution: any constant may be added to one solution to produce another equally
as valid. A variation on this boundary treatment is proposed by Wu et al. in 1D [50] and then also
in 2D [51]. In principle, these approaches are the same as using the Hermite polynomials to extend
the basis, using a bicubically-blended Coons patch in 2D to interpolate the boundary conditions. The
differences are only in terms of the terminology and details of the implementation.
As the developments below illustrate, Sinc methods, like wavelet methods [1], are limited in their direct
application to rectangular or rectangular parallelepiped infinite, semi-infinite or finite geometries. This
is ideal for treating problems with periodic boundary conditions. Attempts to extend Sinc-Galerkin
and similar procedures to more complex geometries either transform the problem as recommended by
McArthur et al. [26] or use a sub-domain (‘building-block’) approach [23]. The former relies on the
conformal invariance of the operator to allow the results to be transformed back to the original prob-
lem [52]. Unfortunately, not every operator of interest satisfies this criterion. The latter divides the
domain into ‘patches’ [23], on which individual Sinc problems are solved. To avoid memory overruns,
alternating algorithms are used, such as the Schwarz-Alternating Method (SAM), which effectively
handle the coupling on the shared boundary. Overlapping methods are also common, where the ap-
proach is similar but with the patches overlapping one another [22, 23], however, additional difficulties
with the formulation of the discrete system means a simple patching approach is generally preferred.
Care is required when using domain decomposition to ensure that the basis functions in each domain
are capable of being matched at the boundary – in lieu of modifying the basis, finite difference approx-
imations were used in [23]. It is worth noting that the inclusion of finite difference approximations at
the boundary has also been used by [33] as an alternative to modifying the basis to allow for Neumann
boundary conditions. However, this is generally seen to degrade the convergence of the method, with
the improved boundary treatment presented by [50] closely resembling the basis extension method
discussed earlier.

2. Green’s Function for the Helmholtz Equation

Time-harmonic acoustic problems, in the absence of a mean flow, may be represented by a Dirichlet
Helmholtz BVP of the form

(∇2 + k 2 )u(x, ω) = f (x, ω) in Ω (1a)

u(x, ω) = g(ω) on ∂S (1b)

This BVP physically describes an acoustic source f radiating sound at a given frequency ω into the
domain where the resulting acoustic fluctuations are governed by the linear wave equation and are
equal to g on the domain boundary. In the most general case we assume the source term to be a
function of spatially compact support

supp [f (x, ω)] = y ∈ ΩS

where ΩS ∈ Ω represents a ‘source region’. One method of developing an analytical solution to this
BVP is to transform the differential equation into an integral equation in which we can substitute
the boundary conditions. These boundary conditions are incorporated in the integral kernel through
the Green’s function Ĝ(x, y, ω). This function represents the outgoing wave produced by an impulsive
unit point source at y subject to the same boundary conditions as the original problem. It is the
solution to the singular form of the original BVP, where our point source is described by the Dirac
delta function Equation (2).
 Green’s Functions using Sinc-Galerkin 5

(∇2 + k 2 )Ĝ(x, y, ω) = δ(x − y) in Ω (2a)

Ĝ(x, y, ω) = g(y, ω) on ∂S (2b)
In the absence of boundary conditions we may consider this wave to propagate in free-space and
hence we often use the free-space Green’s function, readily available in analytical form [13], to solve
the problem. However, the introduction of boundaries and hence associated boundary conditions
complicate matters as an analytical representation of Ĝ is difficult to construct. In the remainder of
this paper we show that the Sinc-Galerkin method can be a robust and accurate means of numerically
computing Green’s functions for homogeneous Dirichlet Helmholtz BVPs within a particular class of
geometries.
We have chosen to solve the 2D Helmholtz equation rather than a higher-dimensional version or indeed
the time-dependent wave equation for simplicity although the Sinc-Galerkin may be easily extended
to higher dimensions with little practical effort. A typical use of the resulting Green’s function would
be to integrate it with an acoustic source term as part of an acoustic analogy. Accurate descriptions
of acoustic sources are usually based on the frequency content as producing a time-dependent source
description from either computational or experimental data or using a mathematical model can be
difficult, especially for real sources which may consist of random behaviour [48]. However, we note that
if we wished to compute a Green’s function for the wave equation rather than the Helmholtz equation,
the Sinc-Galerkin method has the advantage of being able to treat the time variable in the same way
as the space variables; a conformal map is used to map the time functions to the interval [0, ∞]. Hence,
the method may be applied just as easily to the wave equation as the Helmholtz equation without the
need for a complicated time-stepping algorithm or solution of a system of ODEs in time [26]. This
makes the solution of time-dependent problems both quick and stable.

3. 1D Sinc Approximations
We start by presenting the key formulae associated with Sinc numerical methods before applying
extending them to our particular application. The one-dimensional Sinc series for a function f (x) on
the real line x ∈ R may be defined for step size h > 0 as
∞  
X x − kh
C(f, h)(x) ≡ f (kh)sinc (3)
−∞
h
where
(
sin(πz)
πz 6 0
z=
sinc(z) ≡
1 z=0
This series converges for the Paley-Wiener

class of complex functions [43] in which case C is called
the cardinal function of f . Since sinc x−kh
h = 1 when x = kh, the cardinal function can be used
to interpolate f at the collocation points {nh}∞ n=−∞ . Although this cardinal function representation
is exact for Paley-Wiener functions, we may also use the series to accurately approximate functions
outside this class. Furthermore, we truncate Equation (3) in practice such that −M ≤ k ≤ N . In
these circumstances, the 1D sinc interpolation formula for f is stated as [21]

N  
X x − kh
f (x) ≈ f (kh)sinc (4)
h
k=−M

Integration of the series representation Equation (3) over the whole the real line, followed by a
truncation of the resulting series, defines the 1D quadrature rule
6 Harwood and Dupère

Z ∞ N
X
f (x)dx ≈ h f (kh) (5)
−∞ k=−M

When using these approximate formulae, it is important to select the parameters such that the er-
ror in the approximation behaves exponentially. This ensures the series representation remains rapidly
convergent. Accurate approximations may then be achieved using a relatively small number of terms
in the series. This is a particular strength of sinc series representations over similar representations
such as Fourier series. The mathematical rules for parameter selection are established for a number
of problems in [21], although in all cases, knowledge of the solution behaviour at the endpoint of
the interval is required. In practice, such information is rarely available and hence a set of ‘general’
selections are made when performing our calculations which are detailed later.
For many boundary value problems, the domain boundaries are finite, and Equation (4) may be
mapped from x ∈ R to the interval x0 ∈ [a, b] the conformal mapping function φ : x 7→ x0 defined by
 
x−a
φ(x) = ln (6)
b−x

If the inverse mapping function is defined by ψ(x0 ) = φ−1 (x0 ) then we have the new interpolation
formula for the interval [a, b] as

N
X
f (x) ≈ f ◦ ψ(kh) [S(k, h) ◦ φ(x)] (7)
k=−M

as well as the new quadrature formulation

Z b N
X f
f (x)dx ≈ h ◦ ψ(kh) (8)
a φ0
k=−M

where we have introduced a shorthand notation for a translated sinc function defined by
 
z − kh
S(k, h)(z) ≡ sinc (9)
h
Figure 1 shows the basis functions k = −1, 0, 1 on the interval [0, 2] as an example of the translation
and mapping of the basis set. For clarity we have not plotted the curves immediately near the end
points. Had we done so, their rapid decay to zero would be evident.
Finally, as many differential equations contain derivatives we also define the following [15]

(n) dn
δjk ≡ hn [S(j, h) ◦ φ(x)] n = 0, 1, 2, . . . (10)
dφn x=xk

(0) (1)
where xk = φ−1 (kh) = ψ(kh). Specifically for the context of this paper we note the values of δjk , δjk
(2)
and δjk by observing φ(xk )|x=xk = φ φ−1 (kh) ≡ x|x=kh and

dn dn
  
x − jh
[S(j, h) ◦ φ(x)] = sinc
dφn x=xk dxn h x=kh

giving us
 Green’s Functions using Sinc-Galerkin 7

Figure 1. 3 Sinc Basis Functions (shifted to their centring points (circles) and trans-
formed to interval (0,2). The blue curve is k = −1, the green k = 0 and the red k = 1.)

(
(0) 1 j=k
δjk ≡ [S(j, h) ◦ φ(x)]|x=xk = (11a)
0 j 6= k
(
(1) d 0 j=k
δjk ≡ h [S(j, h) ◦ φ(x)] = (−1)k−j (11b)
dφ x=xk k−j j 6= k
( 2
(2) d2
2
− π3 j=k
δjk ≡h [S(j, h) ◦ φ(x)] = −2(−1)k−j (11c)
dφ2 x=xk (k−j)2 j 6= k

Next we introduce the Sinc-Galerkin procedure for a two-dimensional Helmholtz BVP and show the
simple extension of the formulae presented in this section to two dimensions.

4. 2D Sinc-Galerkin for the Helmholtz Equation

The 2D Sinc-Galerkin method for a homogeneous Dirichlet BVP may be derived by considering the
Helmholtz equation on a rectangle (a, b) × (c, d) subject to homogeneous Dirichlet conditions on the
edges Equation (12). Alternatively, as the problem is separable, we could write the original BVP as two
ordinary BVPs: the corresponding BVPs in x for a given yk and in y for a given xj . The advantage
of this latter approach include being able to use the simple discretisation schemes established for
1D operators [3, 16, 26] to build the higher-dimensional solution. Furthermore, a simple change of
variable [21] (Appendix C) may be used to transform the BVP from one with inhomogeneous Dirichlet
boundary conditions to one with homogeneous conditions, as in [42]. As we are considering only the
homogeneous Dirichlet problem, we may treat the problem as a whole.

(∇2 + k 2 )u(x, y) = −f (x, y) (x, y) ∈ Ω = (a, b) × (c, d) (12a)

u(x, y) = 0 (x, y) ∈ ∂S (12b)
Extending Equation (8) to two dimensions we may approximate the solution to Equation (12) as a
sinc series
8 Harwood and Dupère

Nx Ny
X X
uA (x, y) = ujk Sjk (x, y) (13)
j=−Mx k=−My
In this case, the coefficients at the nodal points of each sinc function in the basis set are unknown
since the solution itself is unknown. These unknown are represented by the coefficients ujk . The basis
functions {Sjk (x, y)} use the shorthand of Equation (9) are defined

Sjk (x, y) = [S(j, hx ) ◦ φx (x)] [S(k, hy ) ◦ φy (y)] (14)

The required mapping function for the x-direction is given by
x−a
φx (x) = ln
b−x
which maps (a, b) to (−∞, ∞) and for the y direction
y−c
φy (y) = ln
d−y
which maps (c, d) to (−∞, ∞). We initiate the Galerkin method by considering the orthogonalisation
of the residual R = −(∇2 + k 2 )uA − f with respect to the original basis set of functions Equation (14)

−(∇2 + k 2 )uA − f , Spq = 0

The inner products are defined by

ZZ
hh, gi = h(x, y)g(x, y)w(x)v(y)dxdy
Ω
where w(x)v(y) is a weight function. Green’s second identity
ZZ Z  
2 2

∂h ∂g
g∇ h − h∇ g dxdy = g −h dΓ
Ω ∂S ∂n ∂n
can be used to swap the derivative operations from the solution to the basis functions for which the
derivatives are known (Equation (10)). We therefore get, for a constant wave number, the equivalent
integral equation
ZZ ZZ ZZ
2 2
− uA ∇ gdxdy + BT − k uA gdxdy = f gdxdy
Ω Ω Ω
where BT is the boundary integral. As a means of keeping the notation compact we have used g to
represent the weights and basis functions such that

g ≡ g(x, y)
= gx (x)gy (y)
= [S(p, hx ) ◦ φx (x)] w(x) [S(q, hy ) ◦ φy (y)] v(y)
We continue by expanding the 2D Laplacian operator in the first term in the integral equation and
assume the weight vwqis selected to cause BT to vanish. A suitable choice of weight function to achieve
this is w(x)v(x) = 1/ φ0x (x)φ0y (y). For a more detailed discussion of appropriate weighting functions
see [6]. Our weighted residual statement is, therefore:
ZZ ZZ ZZ ZZ
− uA gy gx00 dxdy − uA gx gy00 dxdy − k 2 uA gdxdy = f gdxdy (15)
Ω Ω Ω Ω
The weighted residual statement is then discretised by substituting the approximation for uA and
applying the 2D quadrature rule
 Green’s Functions using Sinc-Galerkin 9

ZZ Nx Ny
X X z(xj , yk )
z(x, y)dxdy ≈ hx hy
φx (xj )φ0y (yk )
0
j=−Mx k=−My

which is equivalent to Equation (8) extended to two dimensions where z(x, y) is an arbitrary 2D
function and the nodal points are represented by combinations of (xj , yk ). A special case of the
quadrature rule useful for our analysis, for given values of p and q, is [15]
ZZ
G(xp , yq )
G(x, y) [S(p, hx ) ◦ φx (x)] [S(q, hy ) ◦ φy (y)] dxdy ≈ hx hy 0 (16)
φx (xp )φ0y (yq )
(0) (0)
where implicitly the known results of δpj and δqk from Equation (11) have been used. The product
rule is used on the weighted residual statement (Equation (15)) to expand the derivatives in the inner
products. Specifically we note that

00
gx00 (x) = [[S(p, hx ) ◦ φx (x)] w(x)]
0
= [S 0 w + w0 S]
= [S 00 w + S 0 w0 + w00 S + w0 S 0 ]
= [S 00 w + w00 S + 2S 0 w0 ]
but S 0 = S̃ 0 φ0x where S̃ 0 is viewed as a function of φx . Therefore applying the product rule again
2
S 00 = S̃ 0 φ00x + S̃ 00 (φ0x ) gives us
h i
2
= S̃ 00 (φ0x ) w + S̃ 0 (wφ00x + 2w0 φ0x ) + Sw00 (17)

and likewise for gy00 (y). Application of the quadrature rule to the discretised weighted residual statement
gives us the following discretised equation for a given combination of basis functions

Nx
1 (1) φ00x 00
   
v X 1 (2) 0 0 (0) w
(yq ) − 2 δpj φx w − δ w + 2w − δpj 0 (xj )ujq
φ0y hx hx pj φ0x φx
j=−Mx
Ny
1 (1) φ00y 00
   
w X 1 (2) 0 0 (0) v
+ (x p ) − δ φ v − δ v + 2v − δ (yk )upk
φ0x h2y qk y hy qk φ0y qk 0
φy
k=−My
vw f vw
− k2 0 0
(xp , yq )upq = 0 0 (xp , yq ) (18)
φx φy φx φy
Allowing p and q to vary over their range of values we generate a system of equations which after
some slight rearrangement can be expressed as
h iT
D(φ0x )Ã(w)D(φ0x )V + V D(φ0y )Ã(v)D(φ0y ) − k 2 V = G
where D represents a diagonal matrix and
 00
2ω 0 ω 00
  
1 (2) 1 (1) φ
Ã(ω) = − I − I D + − D
h2 h (φ0 )2 φ0 ω (φ0 )2 ω
(19)
V = D(w)UD(v)
G = D(w)FD(v)
where ω can be either v or w and the elements of the matrices I(n) are evaluated using Equation (11).
Using the definitions D(φ0x )Ã(w)D(φ0x ) = Ax and D(φ0y )Ã(v)D(φ0y ) = Ay we represent this system
as a generalised Sylvester equation
10 Harwood and Dupère

Ax V + VATy − k 2 V = G
⇒ (Ax − k 2 Imx )V + VATy = G (20)
Imx is here an identity matrix of size mx × mx .

5. Solution Method
The two most common solution methods for this system are presented in [21] and comprise:
1. Rewrite the system as a large sparse linear system by concatenating the vectors and matrices of
the current system.
2. If we can assume Ax −k 2 Imx and ATy are diagonalisable (depends on the choice of weight function
and order of the problem) we can use factorisation to construct V and hence U.
The first method requires no presumption of the properties of the matrices and is a more tractable
approach when the PDE has non-constant coefficients or contains operators and terms of low order
[26]. Furthermore, using the first method of solution, the unknowns in V are treated as a vector which
lends itself to parallel computing systems. The former method is, therefore, selected over the latter,
although details of the latter are given in Appendix B. One obvious disadvantage of the large sparse
system approach is the memory required to store the large matrices involved. This can put restrictions
on the chosen basis size, and hence approximation accuracy.
Solution using the former method begins by considering the following theorem [15] for concatenating
the individual systems.
Theorem 1. For unknown matrix U, let there be a series of linear systems

A1 UB1 + . . . + An UBn = G
where the dimensions of the matrices Ai and Bi are mx × mx and my × my respectively and therefore
both U and G are mx × my . This is equivalent to

Cco(U) = co(G)
where

C ≡ BT1 ⊗ A1 + . . . + BTn ⊗ A
and the concatenation operator co( ) is defined formally in the appendix of [21].
We may thus rewrite the system in this form with

C = Imy ⊗ Ax − k 2 Imx + Ay ⊗ Imx

Also,

coV ≡ co (D(w)UD(v))
= [D(v) ⊗ D(w)] co(U)
= Dvw co(U)
and hence
co(G) = Dvw co(F)
The Gauss-Jordan elimination method may be applied to increase computational efficiency, exploiting
the special structure of the resulting linear system [25, 26].
 Green’s Functions using Sinc-Galerkin 11

5.1. Parameter Selection

The optimum choice of Mx and Nx depends on the behaviour of the solution at the boundaries of
the interval. Suitable choice enhances the convergence rates of the method as the theorems of [21]
illustrate. However, for practical problems, the solution is likely to be unknown and hence we use the
advice of Mueller and Shores [31] to make sensible selections. From this paper we know
 
α
N= M +1
β
and also
r
πd
h=
αM
and also choose α = β = 1. In addition, we use the conservative choice of d = π/4.

5.2. Removal of the Singularity

At this point it is worth mentioning that the exponential convergence of the Sinc-Galerkin method is
only maintained if the problem has at worst regular singular points at the boundary [43]. By introducing
a singular source within the domain, we violate the conditions for exponential convergence of the sinc
series. However, one way to proceed could be to remove the singularity from the solution and then
apply the Sinc-Galerkin procedure to the non-singular problem, rather than tackling the problem
directly. This, of course, would require knowledge of the nature of the singularity in advance, which
may not be available. If, however, we assume the singularity to be of a logarithmic nature we may
write the solution to Equation (1) as

u(x, y) = ln v(x, y) + ū(x, y)

where ū represents the non-singular part of the solution to be expressed as a sinc expansion. The
equation to solve is then

Lū = −(f + L ln v)
≡ f¯
where L is the Helmholtz operator in the case studied here. Alternatively, we could use classical results
to construct the singularity. For example, the known expansion for a corner singularity is used in Marin
et al. [24]. Such a procedure is not the focus of this paper. Instead we examine how we may directly
apply the Sinc-Galerkin procedure to problems whose solutions contain unknown domain singularities.

6. Numerical Validation
We now apply the 2D Sinc-Galerkin formulation to compute Green’s functions within a rectangular
patch subject to homogeneous Dirichlet-type boundary conditions (g = 0 in Equation (1)). Although
we choose to compute Green’s functions for the Helmholtz equation, the formulation presented here
may also be used for Laplace’s equation without significant modification. As discussed in Section 1, the
use of the Green’s function to solve realistic boundary value problems is a fast and accurate alternative
to common numerical schemes such as FEM or finite differences provided the Green’s function can be
found. Traditionally, the solution of many acoustic BVPs would have been undertaken using an acous-
tic analogy [4, 27, 35]. However, the lack of available Green’s functions for more general geometries
restricts its present use. Choosing the domain to be a patch we allow the present formulation to be
included in a more complicated patching or overlapping algorithm [22, 23]. Such algorithms allow the
method to be applied to a wider class of geometries. Through the selection of appropriate conformal
maps, it is also possible to apply the Sinc-Galerkin method to semi-infinite and infinite rectangular
12 Harwood and Dupère

geometries. In these cases the Sinc-Galerkin method does not require artificial boundaries or non-
reflecting boundary conditions as would be necessary for use of the FEM for the same problem, which
can be particularly useful. We provide the analytical solutions for both the examples below which
we use as a means of assessing the accuracy of our results. Furthermore, the speed of computation
is assessed by comparing the results to those obtained using a constant-element direct BEM imple-
mentation. These examples demonstrate that the Sinc-Galerkin method may be used to numerically
evaluate Green’s functions quickly and with accuracy. The extension of the present formulation to 3D
is straightforward [21] but is not performed here to ensure the results may be clearly understood.
As a means of examining convergence, we vary the mesh density uniformly up to a maximum of 3 times
the lowest mesh density in order to show the trend. We also compare the results to those computed
using a direct-collocation BEM implementation over a range of mesh densities. As before, the ratio
between the largest and smallest mesh density is 3. In addition, we record the execution time in both
cases. All computations are performed using MATLAB on an Intel Core i7 870 CPU.
6.1. Source Specification
In order to configure the Sinc-Galerkin method for our problems (Equation (1)), we must specify
a suitable source function f . As we are considering the BVP defined by Equation (2) along with
homogeneous boundary conditions, analytically, the source function is equal to the Dirac delta func-
tion centred at (y1 , y2 ). A suitable numerical approximation to the delta function would be required
for many numerical methods, usually by considering the well-known limits of a Gaussian function
[2, 30] or a Lorentzian function. Moreover, Tornberg and Bjõrn [46] present an analysis of a variety of
regularisation techniques for a delta function with compact support. However, as we use a Galerkin
weighted residual technique, the application of a numerical quadrature to evaluate the source term
inner product is not strictly necessary. The source term is a delta function which is an integrable
singular function. Therefore we exploit the sifting property to express the source inner product ana-
lytically avoiding any potential error introduced by an approximation. Using the coordinate definition
of the previous section and denoting the source location in the rectangle by (x0 , y0 ) we have
ZZ
hf, Spq i = δ(x − x0 )δ(y − y0 ) [S(p, hx ) ◦ φx (x)] [S(q, hy ) ◦ φy (y)] dxdy
Ω
= Spq vw(x0 , y0 )
In order to incorporate this into Equation (18) we need to manipulate this further; in moving
from Equation (18) to Equation (20) we have divided both sides by hx hy , multiplied through by
φ0x φ0y (xp , yq ) and then, due to the definition of G, multiplied through by vw(xp , yq ). Therefore, we
incorporate the inverse of these operations in the definition of f in Equation (18) giving us a source
function defined by
 
1 vw
fpq (x0 , y0 , xp , yq ) ≡ [Spq vw] (x0 , y0 ) 0 0 (xp , yq )
hx hy φx φy
6.2. Green’s Function for the 2D Helmholtz equation
We consider the Green’s function for the Helmholtz equation in a rectangle subject to homogeneous
Dirichlet boundary conditions on the walls (Equation (2) with g = 0). The analytical solution to the
problem may be constructed using a double Fourier sine transform method thus [38]: Starting with

∇2 + k 2 u(x, y) = −f (x, y)

and letting
f (x, y) = f0 F (x, y)
gives
∇2 + k 2 u(x, y) = −f0 F (x, y)


 Green’s Functions using Sinc-Galerkin 13

Given the boundary
conditions

we may apply double Fourier sine transforms to both sides by multi-
mπy
plying by sin nπx
a sin b and integrating over the rectangle

Za Zb  nπx   mπy  Za Zb  nπx   mπy 

 2
∇ u(x, y) + k 2 u(x, y) sin

sin dydx = −f0 F (x, y) sin sin dydx
a b a b
0 0 0 0
Integrating the left hand side by parts term by term and applying the boundary conditions gives

Za Zb Za Zb
n2 π 2  nπx   mπy  m2 π 2  nπx   mπy 
− 2 u(x, y) sin sin dydx − u(x, y) sin sin dydx
a a b b2 a b
0 0 0 0
Za Zb  nπx   mπy 
2
+k u(x, y) sin sin dydx
a b
0 0
Factorising and rearranging gives

RaRb nπx

mπy


Za Zb  nπx   mπy  f0 F (x, y) sin a sin b dydx
00
u(x, y) sin sin dydx = n2 m2


a b π2 a2 + b2 − k2
0 0
Taking the inverse transform gives for a general source F and (η, ξ) ∈ (0 ≤ x ≤ a, 0 ≤ y ≤ b)

Za Zb 
∞ X ∞  nπx
sin mπy

 
4f0 X sin a b
 nπη  mπξ 
u(x, y) = n2 m2

F (η, ξ) sin sin dξdη
ab n=1 m=1  π 2 2
a2 + b2 − k 0 0
a b 

For the Green’s function we substitute F (η, ξ) = δ(η − η0 )δ(ξ − ξ0 ) and evaluate the integral to give
the following analytical Green’s function for a source located at (η0 , ξ0 )
  
nπx mπy nπη0 mπξ0




∞ ∞
4f0 X X  sin a sin b sin a sin b
u(x, y) = 2 2

 (21)
ab n=1 m=1 π 2 na2 + mb2 − k
2

6.3. Centred Source

The Sinc-Galerkin method is applied to solve the BVP. The number of terms used in the sinc expan-
sion (size of the basis) in both the x and y-directions are the same (as the solution is symmetric) and
equal to 8 uniformly spaced values between M = 12 and M = 36. This range allows us to increase
the mesh density by up to a factor of 3. Similarly, we configure a typical direct BEM implementation
to use a range of 4 to 12 constant elements per unit length (also a factor of 3) to solve the same
problem and compare the results. We also vary the frequency (and hence the wave number k). We test
4 frequencies in the range 20 Hz to 2000 Hz. The lower limit is the lowest audible frequency and the
upper limit chosen so that the plots do not contain too many high order modes making them difficult
to read. The analytical solution Equation (21) is evaluated at each of the four frequencies truncating
the series after 50 terms. The Sinc-Galerkin results for an expansion of 2M + 1 = 65 and the BEM for
44 elements (12 elements per unit length) are both compared to the analytical solution at the tested
frequencies in Figure 2. All plots use the same colour scale.
Figure 2 illustrates that the analytical solution behaves as expected. The characteristic singularity
of the Green’s function at the centre of the point source is visible and the homogeneous Dirichlet
boundary conditions are met. As the frequency is increased the higher harmonics of the allowable
modes appear which is expected to put increased demands on the ability of the numerical schemes to
resolve the solution. By comparing the Sinc-Galerkin solution to the analytical solution, we see that
14 Harwood and Dupère

the singularity is only captured up to a point, with the precision of the Sinc-based approximation in
the immediate vicinity of the source point being lost. This is the principle source of error when using
the Sinc-Galerkin to solve our class of BVP. This may be attributed to the uneven distribution of Sinc
mesh points: our source is located in the centre of the region under consideration. This coincides with
the region most sparsely populated with mesh nodes. Therefore, it is expected that steep gradients
are difficult to represent in this region. In addition, the peak of the sinc function has a value of unity
and hence to capture a singularity requires a large coefficient which would impact the contribution of
this function over the whole of the domain. The resulting behaviour is, therefore, a trade-off between
overall accuracy and the ability to capture the singularity. If the method were included as part of a
domain decomposition approach to the problem, the sub-domain boundaries could be defined such
that the singularity is situated on the boundary. In this case, the exponential convergence of the
method could be restored. Furthermore, although the BEM results indicate that the BEM captures
the singularity more precisely than the Sinc-Galerkin method, there is evidence of a loss of definition
in Figure 2b at the peaks and troughs at the higher frequencies tested. Therefore, qualitatively, the
Sinc-Galerkin captures the general behaviour of the solution with a higher degree of fidelity than the
BEM, the latter demonstrating a notable resolution-driven sensitivity to frequency.
From both the Sinc-Galerkin and BEM results, we also note that the boundary conditions are accu-
rately satisfied by both methods: the BEM allows the boundary conditions to be specified precisely
as in input to the method and the nature of the Sinc function itself means it decays to zero at the
edge of the domain.
We compare the execution time and the relative error in all cases (Figure 4a (log scale) and
Figure 4b). For most of the tested mesh densities tested, the Sinc-Galerkin method performs reasonably
well with the mean relative error remaining less than 15%. The boundary element method sees a
reduction in error with mesh density more steeply than the Sinc-Galerkin method with a 3 fold
increase in the mesh density reducing the error to below 1% for the lowest two frequencies. When
using BEM the largest errors are typically found near the boundary. Although a large relative error
at these points is expected, as the analytical solution tends to zero in these regions, the absolute
error is also larger than the Sinc-Galerkin method in these regions. Due to the nature of internal
point evaluation for the BEM, singular collocation points on the boundary are close by and sensitive
to numerical inaccuracies. In the next section we repeat the analysis of the present section but for
a source located near the left hand boundary of the domain in order to determine how boundary
proximity affects the accuracy of the Sinc-Galerkin method. Consequently, we also observe how this
situation affects the BEM results.
The sensitivity to frequency of the BEM is more dramatic than the Sinc-Galerkin method. The Sinc-
Galerkin method performs poorly at the highest frequency tested initially due to the basis being of
insufficient size to ensure the expansion converges to a suitable degree as to accurately represent the
field. As the number of terms in the expansion is increased, the error begins to reduce. Although at the
highest mesh density tested the relative error is still large, it is driven by a combination of localised
error in the vicinity of the singularity, specifically at the ‘tip’, and small errors in regions where the
analytical representation is near zero resulting in large relative errors. Figure 2 supports this assertion.
In fact, the unusual trend in the relative error is not mirrored in the absolute error results (Figure 3).
In the figure the absolute errors are lower for the BEM for the lowest frequency but relative errors
are comparable. This is largely driven by the errors in the singularity region where the Sinc-Galerkin
fails.
From Figure 4a, the convergence rates for both the BEM and the Sinc-Galerkin are different, with
the BEM converging quicker with a trebling of the mesh density compared with the Sinc-Galerkin.
However, despite the presence of a domain singularity, the lines are roughly linear which suggests an
exponential variation (logarithmic scale). Figure 3 also indicates the limit of exponential convergence
based on the error at the lowest mesh density. This limit is computed from setting the error equal to

√
A(f )e− N
)
 Green’s Functions using Sinc-Galerkin 15

(a) f=20 (left column) and f=680 (right column)

equivalent to setting k = 1 in the expression given in Section 1.1. The constant A(f ) is frequency
dependent and is chosen such that the curves coincide at the lowest mesh density. The convergence is
generally uneven with mesh density and varies with frequency with the lowest frequency converging less
than exponentially with mesh density. However, in the case of the higher frequencies, the convergence
at worst with the exponential convergence limit and in some cases exceeds this. We conclude that the
exponential convergence of the method can still be guaranteed but the decay is significantly slower
than other studies whose solutions do not contain a singularity. This is expected as spectral methods
generally only guarantee rapid convergence when the functions being represented are smooth without
strong singularities or discontinuities. Convergence rates may be improved by implementing the double
exponential transformation ([34, 45, 51]) but is not performed here.
Despite these apparent drawbacks, the Sinc-Galerkin still out performs the BEM in most of the
16 Harwood and Dupère

(b) f=1340 (left column) and f=2000 (right column)

Figure 2. Results for the BEM n=12 (bottom row), analytical expression Equa-
tion (21) (middle row) and Sinc-Galerkin M=32 (top row) for a source centred in the
unit square

combinations tested while still taking less time to execute (Figure 4b). Specifically at lower mesh
resolutions, it can be up to 2.5 times longer to use the BEM over the Sinc-Galerkin. This is due
to the simplicity of the sinc quadrature rule which avoids expensive numerical integration. As the
mesh density increases, however, the time taken scales non-linearly for the Sinc-Galerkin and roughly
linearly for the BEM although with an increased frequency dependence. For very large Sinc-Galerkin
systems, not only does the time taken to compute the Kronecker products increase but the storage
 Green’s Functions using Sinc-Galerkin 17

Figure 3. Absolute Error for a centre source for Sinc-Galerkin (squares) and BEM
(circles) for f=20 (solid blue), f=680 (dashed green), f=1340 (dotted red), f=2000
(dot-dash cyan). The corresponding black lines represent the limit of exponential
convergence. The horizontal axis ticks represent an increase in the mesh density com-
parable in both methods.

also increases in an approximately quadratic sense as the dimensions of the Kronecker matrix is
(mx my × mx my ). In such cases, alternative solution methods may be more efficient (Appendix B).
6.4. Near-Boundary Source
For sources near to the boundary, many boundary integral techniques fail due to the proximity of the
singularity to the collocation points. The steep gradients of the exact solution exacerbates the role
of small errors in integration. With the Sinc-Galerkin method, we effectively locate the singularity
in a region of the domain with a denser sinc node distribution. This distribution is a result of the
non-linear conformal transformation of the evenly-spaced sinc nodes along the real line. We expect to
see an improvement in the Sinc-Galerkin’s ability to capture the singularity but a deterioration in the
BEM’s ability, without special treatment [24], as the internal points approach the singular collocation
nodes on the boundary elements. We repeat the same tests as in the previous section but locate the
source with an x-coordinate corresponding to x/W = 0.1 where W = 1 is the domain width. As before
we plot the surfaces for a Sinc-Galerkin expansion of 2M + 1 = 65, the BEM for 44 elements, and
compare to the analytical solution evaluated using 50 terms of the expansion Equation (21). We also
plot the execution time and the relative error in all cases (Figure 6a and Figure 6b).
The analytical solution plotted in Figure 5 again illustrates the expected behaviour including
the singularity located nearer the boundary. The inclusion of higher-order modes and spreading of
the acoustic wave is visible at f > 680 Hz. Although the Sinc-Galerkin results indicate that the
decay of the solution is well-represented away from the singularity, again the immediate vicinity of the
18 Harwood and Dupère

(a) f=20 (solid blue), f=680 (dashed green), f=1340 (dotted red), f=2000 (dot-dash cyan)

(b) f=20 (solid blue), f=680 (dashed green), f=1340 (dotted red), f=2000 (dot-dash cyan)

Figure 4. Relative error (top) and execution time in seconds (bottom) for Sinc-
Galerkin (squares) and BEM (circles) for the range of mesh densities tested (horizontal
axis) for a centred source.

singularity is not captured in its entirety. Qualitatively, the BEM again captures the singularity well
but struggles at the higher frequencies Figure 5b. The solution is over-predicted at peaks and troughs.
This is more noticeable with the near-boundary source than with the centred source, although this
may be due to the overall solution being an order of magnitude greater in these regions as indicated
by the colour scale. Interestingly, at a frequency of f = 1340 Hz (Figure 5b), the phase of the result
is incorrect by 180 degrees for the BEM result. Further testing over frequencies 1320 ≤ f ≤ 1350
Hz show particularly poor performance by the BEM. Quantitatively, the results are highly inaccurate
throughout the domain. However, it is only at f = 1340 Hz that we get such a noticeable qualitative
error (the phase shift). At this frequency the condition number of the coefficient matrix in the BEM
 Green’s Functions using Sinc-Galerkin 19

(a) f=20 (left column) and f=680 (right column)

linear system is at a local maximum suggesting that the solution of the system will be highly sensitive
to small errors. If we further increase the resolution, this phase shift effect disappears, the condition
number moves away from the maximum, and the solution given by the BEM becomes much more
accurate.
When comparing Figure 4a and Figure 6a, there are few differences to be seen. However, the error
now increases with frequency across the whole range tested unlike for the centred source. There is a
slight increase in overall error in both methods but, as expected, the BEM suffers more significantly
with the error at the highest two frequencies struggling to converge with mesh density.
In contrast, the Sinc-Galerkin converges exponentially again (straight line on the logarithmic graph)
and convergence remains relatively quick regardless of frequency. The exception is the highest fre-
quency tested where the convergence is less predictable. The behaviour is, however, not inexplicable.
20 Harwood and Dupère

(b) f=1340 (left column) and f=2000 (right column)

Figure 5. Results for the BEM n=12 (bottom row), analytical expression Equa-
tion (21) (middle row) and Sinc-Galerkin M=32 (top row) for a source near the
boundary of the unit square

More terms are required for the expansion to represent the rapid changes in the solution at the highest
frequency. Clearly, there is a minimum number of terms required before the cancellation between the
high frequency modes begins to aid accuracy rather than exacerbate inaccuracies. This explains the
‘hump’ in convergence observable in the figure. Despite this behaviour, the robustness of the Sinc-
Galerkin method means that a lack of mesh density manifests itself in a milder way compared with
the BEM.
Figure 6b shows little sensitivity to source location when compared with Figure 4b. This indicates
 Green’s Functions using Sinc-Galerkin 21

(a) f=20 (solid blue), f=680 (dashed green), f=1340 (dotted red), f=2000 (dot-dash cyan)

(b) f=20 (solid blue), f=680 (dashed green), f=1340 (dotted red), f=2000 (dot-dash cyan)

Figure 6. Relative error (top) and execution time in seconds (bottom) for Sinc-
Galerkin (squares) and BEM (circles) for the range of mesh densities tested (horizontal
axis) for a near-boundary source.

that neither of the two source locations selected here are sufficiently close to the boundary as to cause
any computational difficulties which would require more expensive evaluation methods. In summary,
moving the source closer to the boundary has little effect on both the overall error and the convergence
of the Sinc-Galerkin method, unlike the BEM, which makes it significantly more robust.
22 Harwood and Dupère

7. Conclusions
We have applied the 2D Sinc-Galerkin method to numerically evaluate the Green’s function for the
Helmholtz equation in a rectangular patch subject to homogeneous Dirichlet boundary conditions.
The results for a range of frequencies and mesh densities have been compared to a direct BEM
implementation both qualitatively and quantitatively in terms of relative error and execution time.
The results illustrate that increasing the basis size of the Sinc-Galerkin method by a factor of 3
results in a slower convergence than is achieved when increasing the number of boundary elements
by the same factor. Furthermore, unlike the BEM, the domain singularity is not completely captured
in the immediate vicinity of the point itself. However, a suitable means of domain decomposition
could move this singularity onto the boundary and exponential convergence of the method would
be restored. Overall representation of the solution is accurate, robust and quick to evaluate when
compared with the BEM. Only for the largest basis size tested did the time taken to execute storage
requirements exceed those of the BEM. The efficiency of the method may be attributed to the lack
of a requirement for expensive numerical integration. Furthermore, despite the domain singularity,
which reduces the potency of the exponential convergence of the Sinc-Galerkin method, for most
combinations of mesh density and frequency tested, the Sinc-Galerkin achieves lower relative errors
compared to the BEM. Finally, the Sinc-Galerkin is insensitive to a source located near to the boundary
whereas the BEM struggles to maintain accuracy, which highlights the robustness of the method. The
formulation presented may be easily extended to 3D spatial and temporal variations and a formulation
for Neumann (or mixed) boundary conditions is presented.
Acknowledgements
The work was completed as part of the first author’s doctorate degree, supervised by the second author
and funded by the alumni of the University of Manchester through Your Manchester Fund. Some of
this material is to be presented in paper IN14-195 at Internoise 2014 in Melbourne, Australia.

Appendix A. 2D Sinc-Galerkin for Helmholtz Equation with Neumann Boundary

Conditions
For many realistic acoustic problems, boundary conditions involve the spatial derivative of the solution.
This derivative condition is commonly referred to as a Neumann boundary condition. In order to
apply the Sinc-Galerkin method to problems featuring such boundary conditions, we introduce here
a suitable discretisation of the problem.
A.1. Formulation
Two modifications are required to our basis in order to allow their application to the latter BVP.
First, since by definition, 1D sinc functions satisfy homogeneous Dirichlet conditions at the ends of
the interval, they are incapable of representing a non-zero boundary solution. The first modification is
to extend the basis to include functions which represent the solution in the vicinity of the boundary.
Consider the 1D BVP

L [u(x)] = −u00 (x) + p(x)u0 (x) + q(x)u(x) = f (x) x ∈ [a, b]

0
a0 u(a) − a1 u (a) = 0
b0 u(b) + b1 u0 (b) = 0
The values of a0 , a1 , b0 , b1 are positive and if a1 = 0 then b1 6= 0 and vica-versa [21]. The basis may
be extended by two additional functions wa (x) and wb (x) giving approximate solution
N
!
X
uA (x) = u−M −1 wa (x) + uk ξk (x) + uN +1 wb (x)
k=−M
 Green’s Functions using Sinc-Galerkin 23

where ξk (x) are, for now, the usual shifted sinc function basis components Equation (14). The boundary
functions wa and wb are interpolating Hermite functions of order 3 which interpolate the boundary
conditions given in the problem definition as

(x − a)(b − x)2 (2x + b − 3a)(b − x)2

wa (x) = a0 + a1
(b − a)2 (b − a)3
2
(−2x + 3b − a)(x − a) (b − x)(x − a)2
wb (x) = b1 + b0
(b − a)3 (b − a)2

The second modification is to divide the sinc functions in the basis by φ0 to ensure that the
derivative of each basis function is defined at the boundary, which will in turn ensure that the solution
remains defined in these regions. Consider then the Helmholtz BVP with homogeneous Neumann
conditions on a square domain defined as

(∇2 + k 2 )u(x, y) = −f (x, y) (x, y) ∈ Ω = (a, b) × (c, d)

du
=0 (x, y) ∈ Γ
dn

The approximate solution defined in terms of the modified and extended basis is given by

Ny +1 NX
x +1
X
uA (x, y) = ujk ξj (x)ζk (y)
k=−My −1 j=−Mx −1

where

w0 (x)
 j = −Mx − 1
S(j,hx )◦φx (x)
ξj (x) = φ0x (x) ≡ S̃j j = −Mx , . . . , Nx

w1 (x) j = Nx + 1


and

w2 (y)
 k = −My − 1
S(k,hy )◦φy (y)
ζk (y) = φ0y (y) ≡ S̃k k = −My , . . . , Ny

w3 (y) k = Ny + 1


The ensuing method now is based on the following modified inner product orthogonalisation which is
not a true Galerkin procedure but a Petrov-Galerkin procedure since the basis functions used for the
orthogonalisation are a different set to the base used to represent the approximate solution. In fact,
we now orthogonalise the residual with respect to the unmodified basis set

Spq ≡ [S(p, hx ) ◦ φ0x (x)] S(q, hy ) ◦ φ0y (y)

 

for −Mx − 1 ≤ p ≤ Nx + 1 and −My − 1 ≤ q ≤ Ny + 1. We follow the convention [15] and split up
the approximate solution. Our problem produces 9 sub-sets defined thus
24 Harwood and Dupère

Nx
X
uA (x, y) = u(−Mx −1)(−My −1) w0 (x)w2 (y) + uj(−My −1) S̃j (x)w2 (y)
j=−Mx
Ny
X
+ u(Nx +1)(−My −1) w1 (x)w2 (y) + u(−Mx −1)k w0 (x)S̃k (y)
k=−My
Ny Nx Ny
X X X
+ ujk S̃j (x)S̃k (y) + u(Nx +1)k w1 (x)S̃k (y)
k=−My j=−Mx k=−My
Nx
X
+ u(−Mx −1)(Ny +1) w0 (x)w3 (y) + uj(Nt +1) S̃j (x)w3 (y)
j=−Mx
+ u(Nx +1)(Ny +1) w1 (x)w3 (y)

≡ uA1 (x, y) + uA2 (x, y) + uA3 (x, y) + uA4 (x, y) + uA5 (x, y)
+ uA6 (x, y) + uA7 (x, y) + uA8 (x, y) + uA9 (x, y)

We see that uA5 is equivalent to the approximation using for the associated Dirichlet problem (Equa-
tion (13)) with the additional terms arising due to the boundary conditions. We now look at each
inner product of the residual and the unmodified basis set in turn and throughout make use of the
quadrature rule
 (Equation (16)) and the delta notation of Equation (11) and we also note that our
∂2 ∂2
operator L ≡ ∂x2 + ∂y 2 + k 2 . First we start with the inner product where the weighting function
is given by w(x)v(y).

Z Z
hL[uA1 ], Sp (x)Sq (y)i = u(−Mx −1)(−My −1) w000 (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy
ZZ ZZ 
00 2
+ w0 (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy + k w0 (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy

There is no need for any application of Green’s second identity or integration parts, we can just apply
the quadrature rule directly giving us

w000 w

w2 v
hL[uA1 ], Sp (x)Sq (y)i ≈ hx hy (xp ) 0 (yq )u(−Mx −1)(−My −1)
φ0x φy
w00 v

w0 w w0 w w2 v
+ 0 (xp ) 20 (yq )u(−Mx −1)(−My −1) + k 2 0 (xp ) 0 (yq )u(−Mx −1)(−My −1)
φx φy φx φy

Since uA3 , uA7 and uA9 are very similar we can immediately write down their inner product approxi-
mations

w100 w

w2 v
hL[uA3 ], Sp (x)Sq (y)i ≈ hx hy 0
(xp ) 0 (yq )u(Nx +1)(−My −1)
φx φy
w00 v

w1 w w1 w w2 v
+ 0
(xp ) 20 (yq )u(Nx +1)(−My −1) + k 2 0 (xp ) 0 (yq )u(Nx +1)(−My −1)
φx φy φx φy
 Green’s Functions using Sinc-Galerkin 25

w000 w

w3 v
hL[uA7 ], Sp (x)Sq (y)i ≈ hx hy (xp ) 0 (yq )u(−Mx −1)(Ny +1)
φ0x φy
w00 v

w0 w w0 w w3 v
+ 0 (xp ) 30 (yq )u(−Mx −1)(Ny +1) + k 2 0 (xp ) 0 (yq )u(−Mx −1)(Ny +1)
φx φy φx φy

w100 w

w3 v
hL[uA9 ], Sp (x)Sq (y)i ≈ hx hy 0
(xp ) 0 (yq )u(Nx +1)(Ny +1)
φx φy
w00 v

w1 w w1 w w3 v
+ 0
(xp ) 30 (yq )u(Nx +1)(Ny +1) + k 2 0 (xp ) 0 (yq )u(Nx +1)(Ny +1)
φx φy φx φy

Now we consider the inner product associated with uA2 which is

Nx
X Z Z
hL[uA2 ], Sp (x)Sq (y)i = uj(−My −1) S̃j00 (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy
j=−Mx
ZZ ZZ 
+ S̃j (x)w200 (y)Sp (x)Sq (y)w(x)v(y)dxdy + k 2 S̃j (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy

This time, before we apply quadrature we need to apply Green’s second identity to the first inner
product to transfer the derivative onto the sinc basis. Implicitly in the process we assume the weighting
functions will be selected to ensure the boundary terms vanish. The first integral becomes

ZZ ZZ
00
S̃j00 (x)w2 (y)Sp (x)Sq (y)w(x)v(y)dxdy 7→ − S̃j (x)w2 (y)Sq (y)v(y) [Sp (x)w(x)] dxdy

We use the product definition of Equation (17) to arrive at the final approximation.

Nx  00
X w (0) w v w (0) w2 v
hL[uA2 ], Sp (x)Sq (y)i ≈ hx hy uj(−My −1) 0 2
(xj )δpj 20 (yq ) + k 2 0 2 (xj )δpj 0 (yq )
[φx ] φy [φx ] φy
j=−Mx
" (2) (1)
# )
−1 δpj 0 2 δpj 00 0 0 (0) 00 w2 v
[φ ] w + (wφx + 2w φx ) + δpj w (xj ) 0 (yq )
[φ0x ]2 h2x x hx φy

It should be noted that the definition of S̃j ≡ Sj /φ0x has been substituted to ensure the delta definitions
could be used. This accounts for one of the φ0x ’s appearing in the denominator of the first part of the
approximation where the second, which makes the denominator a square, is due to the quadrature
definition. We may immediately write down the expression for uA8

Nx  00
X w (0) w v w (0) w3 v
hL[uA8 ], Sp (x)Sq (y)i ≈ hx hy uj(Ny +1) 0 2
(xj )δpj 30 (yq ) + k 2 0 2 (xj )δpj 0 (yq )
[φx ] φy [φx ] φy
j=−Mx
" (2) (1)
# )
−1 δpj 0 2 δpj 00 0 0 (0) 00 w3 v
[φ ] w + (wφx + 2w φx ) + δpj w (xj ) 0 (yq )
[φ0x ]2 h2x x hx φy

For uA4 and uA6 we have a similar analysis

26 Harwood and Dupère

Ny Z Z
X
hL[uA4 ], Sp (x)Sq (y)i = u(−My −1)k w000 (x)S̃k (y)Sp (x)Sq (y)w(x)v(y)dxdy
k=−My
ZZ ZZ 
+ w0 (x)S̃k00 (y)Sp (x)Sq (y)w(x)v(y)dxdy +k 2
w0 (x)S̃k (y)Sp (x)Sq (y)w(x)v(y)dxdy

As before, we apply Green’s second identity, this time to the second inner product. The second integral
becomes

ZZ ZZ
00
w0 (x)S̃k00 (y)Sp (x)Sq (y)w(x)v(y)dxdy 7→ − w0 (x)S̃k (y)Sp (x)w(x) [Sq (y)v(y)] dxdy

We use the product definition of Equation (17) again to arrive at the final approximation.

Ny
w000 w

X v (0)
hL[uA4 ], Sp (x)Sq (y)i ≈ hx hy u(−Mx −1)k (xp ) 0 2 (yk )δqk
φ0x [φy ]
k=−My
" (2) (1)
#
δqk 0 2 δqk

1 00 0 0

(0) 00 w0 w 2 w0 w v (0)
− 0 2 [φ ] v + vφ + 2v φ + δ v (y k ) (x p ) + k (x p ) (y k )δ
[φy ] h2y y hy y y qk
φ0x φ0x [φ0y ]2 qk

Immediately we find uA6 to be

Ny
w100 w

X v (0)
hL[uA6 ], Sp (x)Sq (y)i ≈ hx hy u(Nx +1)k(xp ) 0 2 (yk )δqk
φ0x [φy ]
k=−My
" (2) (1)
#
δqk 0 2 δqk

1 00 0 0

(0) 00 w1 w 2 w1 w v (0)
− 0 2 [φ ] v + vφy + 2v φy + δqk v (yk ) 0 (xp ) + k (xp ) 0 2 (yk )δqk
[φy ] h2y y hy φx φ0x [φy ]

We then have the formulation for uA5 which is

Ny Nx Z Z
X X
hL[uA5 ], Sp (x)Sq (y)i = ujk S̃j (x)S̃k00 (y)Sp (x)Sq (y)w(x)v(y)dxdy
k=−My j=−Mx
ZZ ZZ 
+ S̃j00 (x)S̃k (y)Sp (x)Sq (y)w(x)v(y)dxdy + k 2 S̃j (x)S̃k (y)Sp (x)Sq (y)w(x)v(y)dxdy

Using Green’s second identity on the first and second integrals this time we get

ZZ
00
S̃j (x)S̃k00 (y)Sp (x)Sq (y)w(x)v(y)dxdy 7→ − S̃j (x)S̃k (y)Sp (x)w(x) [Sq (y)v(y)] dxdy

and

ZZ
00
S̃j00 (x)S̃k (y)Sp (x)Sq (y)w(x)v(y)dxdy 7→ − S̃j (x)S̃k (y)Sq (y)v(y) [Sp (x)w(x)] dxdy

Hence using the quadrature rule and the usual product rule substitution, we get
 Green’s Functions using Sinc-Galerkin 27

Ny Nx 
X X ṽ (0)
hL[uA5 ], Sp (x)Sq (y)i ≈ hx hy ujk (yk )δqk ×
φ0y
k=−My j=−Mx

1 (1) φ00x 00
   
1 (2) (0) w̃
− 2 δpj φ0x w̃ − δpj w̃ + 2 w̃ 0
− δ pj (xj )
hx hx φ0x φ0x
1 (1) φ00y 00
    
w̃ (0) 1 (2) 0 0 (0) ṽ 2 ṽ w̃ (0) (0)
+ 0 (xj )δpj − 2 δqk φy ṽ − δ ṽ + 2ṽ − δqk 0 (yk ) −k 0 0 (xj , yk )δpj δqk
φx hy hy qk φ0y φy φx φy
Inclusion of the usual tilde notation (signifying division by the derivative of the map) as well as
the replacement of all the δ (0) with the appropriate value of coefficient u, is performed to show
the similarity to the Dirichlet formulation (Equation (18)). In particular, note the definitions of the
derivatives
 0
v0 v
ṽ 0 = 0 6=
φ φ0
 00
00 v 00 v
ṽ = 0 6=
φ φ0
Finally, the source function inner product is given by
ZZ
hf (x, y), Sp (x)Sq (y)i = f (x, y)Sp (x)Sq (y)w(x)v(y)dxdy
f wv
≈ hx hy (xp , yq )
φ0x φ0y

A.2. Modified Solution Method

The aim now is to assemble these individual discretised formulations into matrix systems of the form
A1 UB1 + . . . + An UBn . These systems can then be combined to form the generalised Sylvester
equation AUB + CUD = F. If the number of basis functions in the unmodified set with which we are
orthogonalising is given by mx × my then we can say p ∈ mx and q ∈ my . With the original modified
and extended basis, we can identify the basis without the additional extending boundary functions by
nx = mx − 2 and ny = my − 2. Therefore j ∈ nx and k ∈ ny .
First, we consider the sub-systems assembled from the inner products involving uA1 , uA3 , uA7 and uA9
which are all of the form

AUB + ΓU∆ + k 2 ΓUB

where U is a mx × my matrix of coefficients and the remaining matrices are block matrices. The
block matrices in this system are defined in Table 1.
The sub-systems assembled from the inner products involving uA2 and uA8 are all of the form

H̃UB + ΓU∆ + k 2 ΓUB

The block matrices in this system are defined in Table 2.
The sub-systems assembled from the inner products involving uA4 and uA6 are of the form

∆UΓ + BUH̃ + k 2 BUΓ

with the block matrices defined in Table 3.
Finally, the sub-system which forms from the uA5 inner product is of the form

H̃UB + ΓUZ̃ + k 2 ΓUB

28 Harwood and Dupère

Table 1. Block Matrices for Sub-System Type 1

uA1 uA3 uA7 uA9

       
α 0 0 α α 0 0 α
A(mx × mx )
α = w000 w/φ0x (xp ) α = w100 w/φ0x (xp ) α = w000 w/φ0x (xp ) α = w100 w/φ0x (xp )
       
β β 0 0
0 0 β β
B(my × my )
β = w2 v/φ0y (yq ) β = w2 v/φ0y (yq ) β = w3 v/φ0y (yq ) β = w3 v/φ0y (yq )
       
γ 0 γ 0 0 γ 0 γ
Γ(mx × mx )
γ = w0 w/φ0x (xp ) γ = w1 w/φ0x (xp ) γ = w0 w/φ0x (xp ) γ = w1 w/φ0x (xp )
       
δ 0 δ 0
0 δ 0 δ
∆(my × my )
δ = w200 v/φ0y (yq ) δ = w200 v/φ0y (yq ) δ = w300 v/φ0y (yq ) δ = w300 v/φ0y (yq )

Table 2. Block Matrices for Sub-System Type 2

uA2 uA8
   
0 η̃ 0 0 η̃ 0
H̃(mx × mx )
η̃ = HA η̃ = HA

 η  η
H(mx × nx ) (2)
δpj
(1)
δpj
(2)
δpj
(1)
δpj
η= h2x [φ0x ]2 w + hx (wφ00x η= h2x [φ0x ]2 w + hx (wφ00x
 
(0) (0)
+2w0 φ0x ) + δpj w00 (xj ) +2w0 φ0x ) + δpj w00 (xj )
diag[α] diag[α]
A(nx × nx )
α = −1/[φ0x ]2 (xj ) α = −1/[φ0x ]2 (xj )
   
β 0
0 β
B(my × my )
β = w2 v/φ0y (yq ) β = w3 v/φ0y (yq )
   
0 0
 0 diag[γ] 0   0 diag[γ] 0 
0 0
Γ(mx × mx )
γ = w/[φ0x ]2 (xj ) γ = w/[φ0x ]2 (xj )
   
δ 0
0 δ
∆(my × my )
δ = w200 v/φ0y (yq ) δ = w300 v/φ0y (yq )

where the block matrices are defined in Table 4. We may now use Theorem 1 to find the coefficients
contained within the matrix U.
 Green’s Functions using Sinc-Galerkin 29

Table 3. Block Matrices for Sub-System Type 3

uA4 uA6
   
0 0
 η̃   η̃ 
0 0
H̃(my × my )
η̃ = (HA)T η̃ = (HA)T

 η  η
(2) (1) (2) (1)
H(my × ny ) δqk δqk δqk δqk
η= h2y [φ0y ]2 v + hy (vφ00y η= h2y [φ0y ]2 v + hy (vφ00y
 
(0) (0)
+2v 0 φ0y ) + δqk v 00 (yk ) +2v 0 φ0y ) + δqk v 00 (yk )
diag[α] diag[α]
A(ny × ny )
α = −1/[φ0y ]2 (yk ) α = −1/[φ0y ]2 (yk )
   
β 0 0 β
B(mx × mx )
β = w0 w/φ0x (xp ) β = w1 w/φ0x (xp )
   
0 0
 0 diag[γ] 0   0 diag[γ] 0 
0 0
Γ(my × my )
γ = v/[φ0y ]2 (yk ) γ = v/[φ0y ]2 (yk )
   
δ 0 0 δ
∆(mx × mx )
δ = w000 w/φ0x (xp ) δ = w100 w/φ0x (xp )
30 Harwood and Dupère

Table 4. Block Matrices for Sub-System Type 4

uA5
 
H̃(mx × mx ) 0 η̃ 0 η̃ = HA
η
H(mx × nx )  (2) (1)

δpj δpj (0)
η= h2x [φ0x ]2 w + hx (wφ00x + 2w0 φ0x ) + δpj w00 (xj )

A(nx × nx ) diag[α] α = −1/[φ0x ]2 (xj )

 
0
 0 diag[β] 0 
0
B(my × my )
β = v/[φ0y ]2 (yk )
 
0
 0 diag[γ] 0 
0
Γ(mx × mx )
γ = w/[φ0x ]2 (xj )
 
0
 ζ̃ 
0
Z̃(my × my )
ζ̃ = (Z∆)T
ζ
Z(my × ny )  (2) (1)

δqk δqk (0)
ζ= h2y [φ0y ]2 v + hy (vφ00y + 2v 0 φ0y ) + δqk v 00 (yk )

∆(ny × ny ) diag[δ] δ = −1/[φ0y ]2 (yk )

 Green’s Functions using Sinc-Galerkin 31

Appendix B. Alternative Solution Method for the Sylvester Equation

The alternative procedure described in [21, 26] for solving Equation (20) is valid assuming we have
chosen a weight function which enables Ax − k 2 Imx and ATy to be diagonalised. If this is true then
from the classical results there exists matrices P and Q which satisfy

P−1 (Ax − k 2 Imx )P = Λx

Q−1 ATy Q = Λy
where Λx,y is a diagonal matrix of the eigenvalues of Ax − k 2 Imx , ATy and matrices P and Q have
columns corresponding to the associated eigenvectors. This allows us to rewrite the system

Λx W + WΛy = H
where
W = P−1 VQ
H = P−1 GQ
If the spectrum of the matrices are denoted

N
σ(Ax − k 2 Imx ) = {(λx )i }i=−M
x
x
N
σ(Ay ) = {(λy )j }j=−M
y
y

then the components of the unknown matrix W are given by

hij
wij =
(λx )i + (λy )j
where hij refers to the elements of H. Therefore the unknown values of the coefficients stored in U
may be recovered from
   
1 −1 1
U=D PWQ D
v w
where of course we could also have written D v1 ≡ D−1 (v) and D 1
≡ D−1 (w). In order to



w
implement the method, it is necessary to find P, Q, Λx and Λy .

Appendix C. Transformation of Dirichlet Boundary Conditions in 1D

A 1D BVP with inhomogeneous Dirichlet boundary conditions may be transformed into one with
homogeneous boundary conditions by means of the following change of variable, a specific case of the
more general procedure outlined in [21]. Consider the original problem

−v 00 (x) + p(x)v 0 (x) + q(x)v(x) = f (x) x ∈ (a, b)a0 v(a) = a2

b0 v(b) = b2
Using the linear function

l(x) = Ax + B
with coefficients defined by
32 Harwood and Dupère

b0 a2 − a0 b2
A=
D
a0 b2 a − b0 a2 b
B=
D
D = a0 b0 (a − b) 6= 0
Finally u(x) = v(x)−l(x) and f¯(x) = f (x)−L[l(x)] where L is the original differential operator giving

−u00 (x) + p(x)u0 (x) + q(x)u(x) = f¯(x) x ∈ (a, b)

a0 u(a) = 0
b0 u(b) = 0

References
[1] I. Aziz, Islam Siraj ul, and B. Sarler. Wavelets collocation methods for the numerical solution of
elliptic BV problems. Applied Mathematical Modelling, 37(3):676–694, Feb 2013.
[2] Alvin Bayliss and Eli Turkel. Radiation boundary conditions for wave-like equations. Communi-
cations on Pure and Applied Mathematics, 33(6):707–725, 1980.
[3] K. Bowers and J. Lund. Numerical Solution of Singular Poisson Problems via the Sinc-Galerkin
Method. SIAM Journal on Numerical Analysis, 24(1):36–51, 1987.
[4] N. Curle. The Influence of Solid Boundaries upon Aerodynamic Sound. Proceedings of the Royal
Society of London. Series A: Mathematical and Physical Sciences, 231(1187):505–514, 1955.
[5] Mehdi Dehghan and Faezeh Emami-Naeini. The Sinc-collocation and Sinc-Galerkin methods for
solving the two-dimensional Schrödinger equation with nonhomogeneous boundary conditions.
Applied Mathematical Modelling, 37(22):9379–9397, 2013.
[6] M. El-Gamel and A. I. Zayed. Sinc-Galerkin method for solving nonlinear boundary-value prob-
lems. Computers & Mathematics with Applications, 48(9):1285–1298, 2004.
[7] Mohamed El-Gamel. A numerical scheme for solving nonhomogeneous time-dependent problems.
Zeitschrift für angewandte Mathematik und Physik ZAMP, 57(3):369–383, 2006.
[8] Mohamed El-Gamel. A comparison between the Sinc–Galerkin and the modified decomposition
methods for solving two-point boundary-value problems. Journal of Computational Physics,
223(1):369–383, 2007.
[9] Mohamed El-Gamel, S. H. Behiry, and H. Hashish. Numerical method for the solution of spe-
cial nonlinear fourth-order boundary value problems. Applied Mathematics and Computation,
145(2–3):717–734, 2003.
[10] Mohamed El-Gamel and John R. Cannon. On the solution a of second order singularly-perturbed
boundary value problem by the Sinc-Galerkin method. Zeitschrift für angewandte Mathematik
und Physik ZAMP, 56(1):45–58, 2005.
[11] Mohamed El-Gamel, John R. Cannon, and Ahmed I. Zayed. Sinc-Galerkin method for solving
linear sixth-order boundary-value problems. Mathematics of Computation, 73:1325–1343, 2004.
[12] Anna L Hansell, Marta Blangiardo, Lea Fortunato, Sarah Floud, Kees de Hoogh, Daniela Fecht,
Rebecca E Ghosh, Helga E Laszlo, Clare Pearson, Linda Beale, Sean Beevers, John Gulliver,
Nicky Best, Sylvia Richardson, and Paul Elliott. Aircraft noise and cardiovascular disease near
Heathrow airport in London: small area study. BMJ, 347, 2013.
[13] M. S. Howe. Theory of Vortex Sound. Cambridge University Press, Cambridge, 2003.
[14] Z. Y. Huang and W. K. Jiang. An effective method calculating acoustic Green’s function for
closed rectangular cavity using the Ewald’s summation technique. Acta Acustica United with
Acustica, 93(5):853–856, 2007.
[15] Sanoe Koonprasert and Kenneth L. Bowers. The fully Sinc-Galerkin method for time-dependent
boundary conditions. Numerical Methods for Partial Differential Equations, 20(4):494–526, 2004.
 REFERENCES 33

[16] David L. Lewis, John Lund, and Kenneth L. Bowers. The space–time Sinc-Gallerkin method for
parabolic problems. International Journal for Numerical Methods in Engineering, 24(9):1629–
1644, 1987.
[17] M. J. Lighthill. On Sound Generated Aerodynamically. I. General Theory. Proceedings of the
Royal Society of London. Series A: Mathematical and Physical Sciences, 211(1107):564–587, 1952.
[18] M. J. Lighthill. On Sound Generated Aerodynamically. II. Turbulence as a Source of Sound.
Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences,
222(1148):1–32, 1954.
[19] C. Linton. Lattice Sums for the Helmholtz Equation. SIAM Review, 52(4):630–674, 2010.
[20] C. M. Linton. The Green’s Function for the Two-Dimensional Helmholtz Equation in Periodic
Domains. Journal of Engineering Mathematics, 33(4):377–401, 1998.
[21] John Lund and Kenneth L. Bowers. Sinc Methods for Quadrature and Differential Equations.
Society for Industrial and Applied Mathematics, Philadelphia, 1992.
[22] Nancy J. Lybeck and Kenneth L. Bowers. The Sinc-Galerkin Schwarz Alternating Method for
Poisson’s Equation. In Computation and Control IV, volume 20 of Progress in Systems and
Control Theory, pages 247–256. Birkhäuser Boston, 1995.
[23] Nancy J. Lybeck and Kenneth L. Bowers. Domain decomposition in conjunction with sinc meth-
ods for Poisson’s equation. Numerical Methods for Partial Differential Equations, 12(4):461–487,
1996.
[24] L. Marin, D. Lesnic, and V. Mantič. Treatment of singularities in Helmholtz-type equations using
the boundary element method. Journal of Sound and Vibration, 278(1–2):39–62, 2004.
[25] Kelly M. McArthur, Kenneth L. Bowers, and John Lund. Numerical implementation of the
Sinc-Galerkin method for second-order hyperbolic equations. Numerical Methods for Partial
Differential Equations, 3(3):169–185, 1987.
[26] Kelly M. McArthur, Kenneth L. Bowers, and John R. Lund. The Sinc method in multiple space
dimensions: Model problems. Numerische Mathematik, 56(8):789–816, 1990.
[27] W. Möhring. On vortex sound at low Mach number. Journal of Fluid Mechanics, 85(04):685–691,
1978.
[28] Adel Mohsen and Mohamed El-Gamel. On the Galerkin and collocation methods for two-
point boundary value problems using sinc bases. Computers & Mathematics with Applications,
56(4):930–941, 2008.
[29] A. Moroz. Quasi-periodic Green’s functions of the Helmholtz and Laplace equations. Journal of
Physics A: Mathematical and General, 39:11247, 2006.
[30] Philip J. Morris. The scattering of sound from a spatially distributed axisymmetric cylindrical
source by a circular cylinder. The Journal of the Acoustical Society of America, 97(5):2651–2656,
1995.
[31] J. L. Mueller and T. S. Shores. A new sinc-galerkin method for convection-diffusion equations
with mixed boundary conditions. Computers & Mathematics with Applications, 47(4–5):803–822,
2004.
[32] Susheela Narasimhan, Kuan Chen, and Frank Stenger. A Harmonic-Sinc Solution of the Laplace
Equation for Problems with Singularities and Semi-Infinite Domains. Numerical Heat Transfer,
Part B: Fundamentals, 33(4):433–450, 1998. doi: 10.1080/10407799808915042.
[33] Susheela Narasimhan, Joseph Majdalani, and Frank Stenger. A First Step in Applying the Sinc
Collocation Method to the non-linear Navier-Stokes Equations. Numerical Heat Transfer, Part
B: Fundamentals, 41(5):447–462, 2002. doi: 10.1080/104077902753725902.
[34] Ahniyaz Nurmuhammad, Mayinur Muhammad, and Masatake Mori. Sinc–Galerkin method based
on the DE transformation for the boundary value problem of fourth-order ODE. Journal of
Computational and Applied Mathematics, 206(1):17–26, 2007.
[35] F. Obermeier. The influence of solid bodies on low Mach number vortex sound. Journal of Sound
and Vibration, 72(1):39–49, 1980.
34 REFERENCES

[36] Aydin Secer and Muhammet Kurulay. The sinc-Galerkin method and its applications on singular
Dirichlet-type boundary value problems. Boundary Value Problems, 2012(1):126, 2012.
[37] Aydin Secer, Muhammet Kurulay, Mustafa Bayram, and Mehmet Akinlar. An efficient computer
application of the sinc-Galerkin approximation for nonlinear boundary value problems. Boundary
Value Problems, 2012(1):117, 2012.
[38] A.P.S. Selvadurai. Partial Differential Equations in Mechanics: Vol. 1 Fundamentals, Laplace’s
equation, diffusion equation, wave equation., volume 1. Springer-Verlag, New York, 2000.
[39] B. A. Singer, D. P. Lockard, and G. M. Lilley. Hybrid acoustic predictions. Computers &
Mathematics with Applications, 46(4):647–669, 2003.
[40] Wesley C. H. Slemp and Rakesh K. Kapania. Imposing boundary conditions in Sinc method
using highest derivative approximation. Journal of Computational and Applied Mathematics,
230(2):371–392, 2009.
[41] R. Smith, G. Bogar, K. Bowers, and J. Lund. The Sinc-Galerkin Method for Fourth-Order
Differential Equations. SIAM Journal on Numerical Analysis, 28(3):760–788, 1991.
[42] Ralph C. Smith, Kenneth L. Bowers, and John Lund. A fully Sinc-Galerkin method for Eu-
ler–Bernoulli beam models. Numerical Methods for Partial Differential Equations, 8(2):171–202,
1992.
[43] Frank Stenger. A ’Sinc-Galerkin’ method of solution of boundary value problems. Mathematics
of Computation, 33(145):85–109, 1979.
[44] Frank Stenger. Summary of sinc numerical methods. Journal of Computational and Applied
Mathematics, 121(1-2):379–420, 2000.
[45] Masaaki Sugihara and Takayasu Matsuo. Recent developments of the Sinc numerical methods.
Journal of Computational and Applied Mathematics, 164 & 165(0):673–689, 2004.
[46] Anna-Karin Tornberg and Björn Engquist. Numerical approximations of singular source terms
in differential equations. Journal of Computational Physics, 200(2):462–488, 11/1/ 2004.
[47] Meng Wang, Jonathan B. Freund, and Sanjiva K. Lele. Computational prediction of flow-
generated sound. Annual Review of Fluid Mechanics, 38(1):483–512, 2006.
[48] J. Whitmire and S. Sarkar. The computation of flow-generated sound using an acoustic analogy.
AIAA Paper (95-038), 1995.
[49] E. T. Whittaker. On the Functions Which are Represented by the Expansions of the Interpolation
Theory. Proceedings of the Royal Society of Edinburgh, 35:181–194, 1915.
[50] Xionghua Wu, Wenbin Kong, and Chen Li. Sinc collocation method with boundary treatment
for two-point boundary value problems. Journal of Computational and Applied Mathematics,
196(1):229–240, 2006.
[51] Xionghua Wu, Chen Li, and Wenbin Kong. A Sinc-collocation method with boundary treatment
for two-dimensional elliptic boundary value problems. Journal of Computational and Applied
Mathematics, 196(1):58–69, 2006.
[52] Toshihiro Yamamoto. Toward the sinc-Galerkin method for the Poisson problem in one type of
curvilinear coordinate domain. Electronic Transactions on Numerical Analysis, 23:63–75, 2006.
[53] M. Zarebnia and M. Sajjadian. The sinc–Galerkin method for solving Troesch’s problem. Math-
ematical and Computer Modelling, 56(9–10):218–228, 2012.

Adrian R. G. Harwood*
School of Mechanical, Aerospace & Civil Engineering
The University of Manchester
Sackville Street
Manchester
M1 3BB, United Kingdom
e-mail: [email protected]
 REFERENCES 35

Iain D. J. Dupère
School of Mechanical, Aerospace & Civil Engineering
The University of Manchester
Sackville Street
Manchester
M1 3BB, United Kingdom
e-mail: [email protected]