Boundary Element Method or Boundary Integral Method Computational Aspects I have to solve a Helmholtz equation inside a simply connected domain. I know that in general the boundary integral can be written as,
$$\phi(x)=\int_V G(x,x') \rho(x')\ d^3x'+\int_S \left[\phi(x')\nabla' G(x,x')-G(x,x')\nabla'\phi(x')\right] \cdot d\hat\sigma'$$
Since my problem is to solve the Helmholtz equation with the Dirichlet boundary condition the boundary integral can be written as ($\rho$ is zero inside the boundary and $\phi(x)=0$ on the boundary), 
$$\phi(x)=-\int_SG(x,x') \nabla'\phi(x')  \cdot d\hat\sigma'$$
My problem is, how to implement this boundary integral numerically in order to get the eigenvalues and eigenfunctions inside the entire boundary. How can I transform this into a matrix form to obtain the eigenvalues and eigenfunctions?
By the way, the domain is two-dimensional and I am solving this equation in connection with quantum chaos (Quantum chaotic billiard). 
 A: I think the most powerful approach to the general problem of arbitrary domains is the approach detailed in Reviving the Method of Particular Solutions by Timo Betcke and Lloyd N. Trefethen from 2005. [doi] [pdf].  In it they describe a modern modification of the historical method of numerically finding the solution to the helmhotz equation from Approximations and Bounds for Eigenvalues of Elliptic Operators by Fox, Henrici and Moler from 1967. [doi] (fun fact: Moler would go on to create MATLAB)
For completeness, I'll summarize the method.  Philosophically, it is very similar to the variational approaches used to find bounds on energies of eigenfunctions in quantum mechanics.  The basic idea is to choose a possible basis for your eigenfunctions, and then just find a linear combination of those basis functions that satisfies the boundary conditions.  That is, assume we have some basis of parametrized functions:
$$ \phi_{k} (x; \lambda) $$
we will look for solutions to the Helmholtz problem of the form:
$$ \psi(x) = \sum_k c_k \phi_{k}(x; \lambda) $$
and we do this by choosing a set of points on the boundary, for which we will ensure that our solution vanishes, and a set of points on the interior for which the solution doesn't vanish.
The details are worked out in the paper, but this means you choose a set of points $x_i$ on both the boundary and interior, then we define $\sigma(\lambda)$ in the following way:  


*

*Form the matrix
$$ A_{ik}(\lambda) = \phi_k(x_i; \lambda) $$

*takes its QR decomposition
$$ Q_{il}R_{lk} = A_{ik}(\lambda) $$

*then take the smallest singular value of the part of $Q$ for those points on the boundary.


Then, to find your eigenvalues, simply look for minima of $\sigma(\lambda)$ using your favorite numeric scalar minimizer.  Once you've found the eigenvalues, you can find the coefficients determining your solution as the solution to
$$ R c = \hat v $$
where $R$ is the $R$ from our QR decomposition above, $c$ is the vector of coefficients determining your solution and $\hat v$ is the right singular vector of the part of $Q$ defined on the boundary associated with the smallest singular value.
As a form of demonstration, we will consider the classic $L$ shaped domain, the solution to which was first given in the 1967 paper, and which serves as the logo for MATLAB.
To make our method perform well, we will define our basis as:
$$ \phi_k(x; \lambda) = j_{\alpha k} ( \sqrt{\lambda} r ) \sin(\alpha \theta k ) $$
defined from the interior corner of the $L$.  These are the solutions to the helmholtz equation on a semiinfinite wedge, and so they are guarenteed to satisfy our boundary conditions on the interior walls, which will help us out.  So, for our trial points, I'll just take a bunch of random points on the boundary and interior (though I won't put any on the interior walls as I know those are already taken care of by my parametrization)

And then, if we calculate $\sigma(\lambda)$ we get:

Which shows us where we expect to find the eigenvalues, using a scalar minimizer, and finding the eigenfunctions as described, we find for the lowest 6 modes:

Where, for instance our first eigenvalue comes out accurate to 11 decimal places, using only 72 points on the interior and boundary each, and only taking ~ 3 seconds to find all 6 of the lowest energy eigenvalues and eigenfunctions on my machine.  Pretty good if you ask me.  For what it's worth, I've also made the code I used to generate all of these figures available as an ipython notebook where you can see the method implemented.
A: May I understand that, since the domain is 2D, the surface integral is actually the following: $\vec{E}\cdot d\vec{\sigma}' = \vec{E}\cdot\hat{n}dl$, where $\vec{E}$ is an arbitrary 2D vector, $\hat{n}$ denotes the normal of the domain boundary and $dl$ is the line element ? Suppose I understand you correctly and then I propose to switch to the momentum space. 
Let $\phi(\vec{x}) = \sum_k e^{i\vec{k}\vec{x}}\phi_{\vec{k}}$. On substitution , we get 
$$\phi_\vec{k} = \sum_{\vec{k}'} F(\vec{k},\vec{k'})~\phi_{\vec{k}'},$$ where the matrix $F$ is given by $$F(\vec{k},\vec{k}') = -\frac{i}{N}\int_{body}d^2\vec{x}\int_{boundary} dl' e^{i(\vec{k}\vec{x}-\vec{k}'\vec{x}')}\cdot(\vec{k}'\hat{n}')\cdot G(\vec{x},\vec{x}').$$ $N$ denotes the number of points in k-space. I think this formula allows numerical computation.
