1 Introduction

In this paper a computational method for solving the interior Helmholtz eigenvalue problem is explored. The problem is that of finding the values of the wavenumber k and a non-trivial scalar function j such that the Helmholtz equation

Ñ² j(p) + k² j(p) = 0 (p Î D)

(1)

is satisfied in an interior domain D with boundary S and subject to a homogeneous boundary condition of the form

a(p) j(p) + b(p)

¶j(p)

¶n_p

= 0 (p Î S)

(2)

where a(p) and b(p) are known complex-valued functions of p ( Î S) and n_p is the unit outward normal to the boundary at p. The non-trivial solutions k=k^* and j(p) = j^*(p) (p Î D ÈS) are termed the eigenfrequencies and eigenfunctions and they are dependent on the boundary S and the boundary functions a(p) and b(p) and the eigenfrequencies are all real numbers. In this paper we mainly consider the Dirichlet eigenproblem (a(p)=1 and b(p)=0) and the Neumann eigenproblem (a(p)=0 and b(p)=1). Kuttler and Sigillito [1] is a useful background reference for the Helmholtz eigenvalue problem.

A technique for solving this problem has several applications. For example, in acoustics when an eigenfrequency analysis of an enclosure containing an acoustic fluid is required. Such a technique would also be useful when attempting to solve exterior Helmholtz problems via boundary element methods derived in the standard way from the Helmholtz-Kirchoff integral equation or an integral equation arising from expressing the solution as a single- or double-layer potential. For these equations have either no solution or no unique solution at eigenfrequencies of a corresponding interior problem [2], [3]. Thus it is found that methods derived in the standard way from these equations perform poorly in the neighbourhood of these eigenfrequencies [4], [5]. In such circumstances, a technique for locating the problem frequencies would clearly be useful.

The integral operators that arise in corresponding exterior and interior problems are similar. Hence much of the analysis of the integral operators of the exterior problem is applicable to the interior problem and vice-versa. References [6], [7] [8], and [9] consider the condition number of the Helmholtz integral operators as a function of real k. These provide useful background references for the work in this paper as the peaks in condition number occur at eigenfrequencies of the interior Helmholtz equation.

The Helmholtz eigenvalue problem is amenable to solution via finite element or finite difference methods. In these cases, the problem reduces to that of solving a generalised linear eigenvalue problem of the form

(K - k² M) x = 0

(3)

where the matrices K and M (termed the stiffness and mass matrices) in (3) are sparse and structured and are independent of k. Standard computational algorithms are available for solving generalised linear eigenvalue problems. Indeed special techniques (such as iterative methods) are available for solving (3), given the special structure of the matrices and the fact that only a fraction of the full set of eigenvalues are generally required [10]. Hence eigenfrequency analysis of the Helmholtz problem via the finite element or finite difference method is straightforward.

In cases where it is applicable, it is well known that the boundary element method has an important advantage over the finite element and finite difference methods: the partial differential equation governing the domain is reduced to an integral equation relating values of j and [(¶j)/(¶n)] on the boundary only. Hence the dimension of the problem is effectively reduced by one. However, the application of the boundary element method reduces the Helmholtz eigenvalue problem to that of solving an eigenproblem of the form

A_k m = 0

(4)

where the matrix A_k is generally full, having no particular structure but with each component being a continuously differentiable complex-valued function of k.

Because of the main advantage of the boundary element method over finite element and finite difference methods stated earlier, the matrix in (4) is generally much smaller than the matrices in (3), for any given Helmholtz problem and a given level of required accuracy. The disadvantages of this approach are that the eigenvalue problem (4) is non-linear and the components of the A_k matrix are defined in terms of integrals and hence may be costly to evaluate. The solution of non-linear eigenvalue problems are considered in references [11], [12] and [13]. Unfortunately, standard algorithms for solving non-linear eigenvalue problems are not generally available. Hence the application of the boundary element method to the Helmholtz eigenvalue problem is not straightforward.

The problem of solving the Helmholtz eigenvalue problem via boundary element-type methods have been given some consideration by researchers. For example iterative methods such as the secant method are applied to the problem of finding the roots of the equation det (A_k ) = 0 in references [14], [15], [16], [17] and [18]. However, this is not a satisfactory method when the matrix A_k is large [13]. A similar method, based on finding the values of k for which the smallest eigenvalues of A_k is zero is considered in [19]. Unfortunately, these methods are unwieldy since they do not compute the solutions simultaneously, they require a starting point to be chosen for each required eigenfrequency.

In reference [20] a hybrid of the boundary element and finite element method is introduced. The method seems to have the advantage of the finite element method in that a linear eigenvalue problem results and the advantage of the boundary element method in that a solution on the boundary only is sought in the main computation. The method is considered further in references [21], [22].

In general, both eigenfrequencies and eigenfunctions of the Helmholtz problem will be of interest, though in this paper we are principally concerned with results from the computation of the eigenfrequencies. The method considered in this paper is that of approximating each component of the matrix A_k by a polynomial in k in some given sub-range of the full wavenumber range. This allows us to re-write the non-linear eigenvalue problem (4) in the form of a standard generalised eigenvalue problem. Thus all of the eigenvalues in the sub-range are computed simultaneously. The method is applied to the axisymmetric three-dimensional problem where the surface is a sphere and a two-dimensional problem where the boundary is a square. The effectiveness of the method is studied through considering the results of varying the width of the subrange of k, the number of boundary elements and the degree of the interpolating polynomial. (2)