# 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).

• You're interested in how to solve the problem for an arbitrarily shaped domain? Aug 5 '14 at 18:00

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:

1. Form the matrix $$A_{ik}(\lambda) = \phi_k(x_i; \lambda)$$
2. takes its QR decomposition $$Q_{il}R_{lk} = A_{ik}(\lambda)$$
3. 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.

• Dear alemi thank you very much for your detailed answer. This is really helpful. You even have written a code for this question !!! Thanks a lot man. Aug 6 '14 at 14:09
• @SijoJoseph My pleasure. I've been meaning to really understand this method for a while now, this was a great motivation to finally try it out for myself. Also, just so you know, unless you click the little "+50" button next to an answer explicitly, the bounty wont necessarily be awarded automatically. ;) Aug 6 '14 at 14:36
• got it ..I have done it ! Aug 6 '14 at 14:43

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.

• Ok thank you for the reply. How will I determine the eigenfunctions and eigenvalues from this step ? Jul 31 '14 at 19:33
• So, you first construct the F matrix by knowledge of the Green's function G(x,x'). Then you solve the equation for $\phi_{\vec{k}}$, from which you find $\phi(\vec{x})$. Are these the eigenfunctions you are looking for?
– hyd
Aug 1 '14 at 1:16
• Thank you for the answer. It seems that it will work.I accepted it as an answer.But I need to learn it more. By the way could you suggest me some books/articles which deals with this method ? Aug 1 '14 at 16:06
• Sorry I'm not aware of books/papers on this. The idea came out of my experience with Fourier transforms.
– hyd
Aug 2 '14 at 11:59
• No problem. I have figured out the method. I have to compute the Fredholm determinant of a matrix by discretizing the integral equation. That will give me the eigenvalues, from these eigenvalues I can compute the eigenfunction. Anyways thank you for the answer. I was not able to find this method based on the fourrier transform. Aug 2 '14 at 17:12