5  Solution of the Non-linear Eigenvalue Problem

The discrete eigenvalue problems are each of the form (4). The method employed for solving (4) requires that in an interval [kA, kB] of values of k the matrix Ak is approximated by a matrix polynomial in k
Ak » A[0] + k A[1] + ... + km A[m]   for k real.
(9)
Similar methods of approximation are considered in references [24], [25] and [26], but with regard to reducing the computational cost in the standard boundary element solution of Helmholtz problems over a range of wavenumbers. These references do not consider the eigenvalue problem.
The non-linear eigenvalue problem (4) can be replaced with the following eigenvalue problem
[A[0] + k A[1] + ... + km A[m]]m = 0.
(10)
The solutions of (10) are the same as those of the following generalised linear eigenvalue problem

é
ê
ê
ê
ê
ê
ê
ë
A[0]
A[1]
A[2]
.
.
A[m-2]
A[m-1]
0
I
0
.
.
0
0
:
:
:
:
:
:
:
0
0
0
.
.
I
0
0
0
0
.
.
0
I
ù
ú
ú
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ê
ê
ë
m
k m
:
km-2 m
km-1 m
ù
ú
ú
ú
ú
ú
ú
û
(11)

               =   k  é
ê
ê
ê
ê
ê
ê
ë
0
0
0
.
.
0
-A[m]
I
0
0
.
.
0
0
:
:
:
:
:
:
:
0
0
0
.
.
0
0
0
0
0
.
.
I
0
ù
ú
ú
ú
ú
ú
ú
û
é
ê
ê
ê
ê
ê
ê
ë
m
k m
:
km-2 m
km-1 m
ù
ú
ú
ú
ú
ú
ú
û
  .
Equation (11) is amenable to solution by the QZ algorithm [27], which may be invoked, for example, via NAG routine F02GJF [28]. Methods for solving problems of the form (10) are considered in references [29], [11], [30] and [13]. The authors are not aware of any special methods of solving (11) which takes advantage of its special structure or the sparsity of the matrices.
Since the eigenvalues k* of the underlying Helmholtz problem are all real and we are interested only in positive values, it is sufficient to compute the interpolant (9) for the positive real numbers k. However, due to the application of the boundary element method and the approximation (9), the solutions of (11) that correspond to the solution of the underlying Helmholtz eigenvalue problem will tend to have small imaginary parts.
The generalised eigenvalue problem (11) will generally have m ×n solutions. Half of these can be immediately discounted since the eigenvalues occur in pairs [^k], -[^k]. The full set of solutions will contain approximations to the true eigenvalues of the underlying Helmholtz problem. However, many spurious solutions are generally produced as a result of the collocation method and approximation (9). These spurious eigenvalues do not have small imaginary parts and hence they can be sorted from the true eigenvalues. Approximations to the true eigenvalues lying outside the range [kA, kB] may also be produced. These approximations will generally be poor and they can be excluded from the results.
Let [^k], [^(m)] be a typical non-spurious solution to (11). The eigenvalue [^k] is an approximations to the eigenfrequencies k* of the Helmholtz problem. The approximation to the eigenfunctions in D ÈS can be recovered through the substitution of the approximation [^(m)] for s in equation (5) or [^(m)] for j (Dirichlet) or [(j)/(n)] (Neumann) in equations (7)-(8).
In this paper the method employed for deriving the polynomial approximation (9) involves computing Ak at the m+1 Chebyshev (¥ norm) interpolation points for any selected range [kA, kB]. The coefficient matrices A[0], A[1], ..., A[m] in (9) are obtained through Newton's divided differences using the value of Ak at the selected values of k in [kA, kB]. The generalised eigenvalue problem (11) is then solved through invoking NAG routine F02GJF.