4  Boundary Element Method

The boundary element form of the Helmholtz equation is derived through recasting the partial differential equation as a boundary integral equation and applying an integral equation method. The integral equation that is selected for the purposes of this work is the one derived through Green's second theorem. Collocation is applied to derive its discrete form.

4.1  Formulation of the Helmholtz Eigenvalue Problem

The application of Green's second theorem to the Helmholtz equation gives the following equations:
{ Mk j}S (p) + j(p) = { Lk j

n
}S (p)    (p Î D) ,
(14)

{ Mk j}S (p) + c(p) j(p) = { Lk j

n
}S (p)     (p Î S) ,
(15)
where the Helmholtz integral operators, Lk and Mk, are defined as follows:
{ Lk m}P(p) º
ó
õ
P 
 Gk(p,qm(q)  dSq     (p Î D ÈS ),
(16)

{ Mk m}P(p) º
ó
õ
P 
  Gk

nq
(p,qm(q) dSq     (p Î D ÈS ),
(17)
where nq is a unit outward normal to the boundary S at q and p, P Í S , and m(q) is a bounded function defined for q Î P. Gk(p,q) is the free-space Green's function for the Helmholtz equation:
Gk(p,q) = eikr

4 pr
     in three dimensions,
where r=p - q, r=|r|. The function c(p) is defined so that 4 pc(p)  (p Î S) is the solid angle subtended by the exterior at p. Note that if S is smooth at p then c(p) = [1/2].
In order to derive the integral equation form of the eigenvalue problem, the homogeneous boundary condition is enforced on the integral equation formulations (14), (15). The ensuing formulation for the Neumann eigenproblem is as follows:
       { Mk j}S (p) + c(p) j(p) = 0    (p Î S) ,
(18)

       j(p) = - { Mk j}S (p)     (p Î D)  .
(19)
The strategy for solving the eigenvalue problem is thus that of computing the solutions k* and j*(p)  (p Î S) to the integral equation (18). The substitution of j*(p)  (p Î S) into (19) will then yield the domain solution j*(p)  (p Î D).

4.2  Application of Collocation

Collocation is applied to derive the boundary element form of the eigenvalue problems. Let z1(p), z2(p), ..., zn(p) (p Î S) be a linearly independent set of basis functions with the property
n
å
j=1 
zj(p) º 1     (p Î S).
(20)
Let p1, p2,..., pn be the collocation points so that
zi(pj) = dij    (i,j = 1,2,...,n) .
(21)
For the methods considered in this paper, the boundary is smooth at the collocation points, so that c(pj) = [1/2]. This reduces the integral equation formulation (18) of the eigenvalue problems to the following
Ak j » 0 ,
where the matrix Ak is defined as follows
[Ak]ij = {Mk zi }S (pj) + 1

2
dij .
(22)
Approximations to the solutions can now be found through the solution of the following non-linear eigenvalue problem
Ak 
^
j
 
= 0 .
(23)

4.3  Solution of the Non-linear Eigenvalue Problem

For a given problem and a given level of required accuracy, the system of equations (23) will be much smaller than the system of equations arising from the finite element method (13). On the other hand, the straightforward application of the boundary element method results in a non-linear eigenvalue problem (23) and the components of the matrix Ak are defined in terms of integrals and hence may be costly to evaluate.
A typical way of computing the approximations to the eigenfrequencies is to compute the zeros of |Ak|, as considered in [9-12,14]. In Ya Yan Lu and Shing-Tung Yau [27] solutions are sought through finding the values of k at which Ak Ak* has a zero eigenvalue (the * denotes the Hermitian transformation) and in Bai [28] the values of k are sought for which Ak has a zero singular value. Each of these method is likely to be costly since several matrix evaluations (each component computed through the evaluation of an integral) will be necessary (although the matrix components could be approximated by a simple formula such as a polynomial in k, which is cheaper to evaluate). A further drawback of these methods is that the solutions are not computed simultaneously, a starting point is required for each eigenfrequency. Finally, the methods may not readily yield the eigenfunction.
The matrix Ak in (23) has the useful property that its components are continuously differentiable complex-valued functions of the wavenumber k. This means that the components can be satisfactorily approximated by a piecewise polynomial. The approximation of Ak by a matrix polynomial allows us to replace the eigenvalue problem (23) with a standard generalised eigenvalue problem. This method will be described fully in section 6.