The Boundary Element Method

The boundary element method (BEM) is a method for solving differential equations, especially partial differential equations, by discretizing only the boundary of the region. In two dimensions the boundary is a curve; in one dimension the boundary is just the two endpoints of an interval.

Most other numerical methods (finite differences FDM, finite elements FEM, finite volumes FVM, …) discretize the whole region: a grid or a mesh of elements covers the entire domain. The boundary element method discretizes only the boundary, so the resulting linear system is much smaller. This is the main advantage of the method.

The steps of the method are:

1. Transform the PDE in $\mathrm{\Omega }$ into an integral equation on the boundary $\partial \mathrm{\Omega }$.
2. Discretize the boundary $\partial \mathrm{\Omega }$.
3. Solve the integral equation on $\partial \mathrm{\Omega }$ by discretizing it.
4. The result is a linear system $A\mathrm{𝐮}=\mathrm{𝐠}$, where $\mathrm{𝐮}$ is the vector of unknown boundary values.
5. Use the boundary values to obtain the solution at an interior point.

The tool behind the first step is the weighted residual method (WRM), a general way of finding approximate solutions of differential equations. The next sections develop it on simple examples.

---

Weighted Residual Methods

The residual

We want the residual to be as small as possible. Several methods make this precise. All of them start with an approximate form for the solution containing unknown coefficients. Here we take

$$\hat{y}(x)=a{x}^2+b{x}^3.$$

There is no constant term because $y(0)=0$; there is no linear term either, so every choice of $a,b$ satisfies the initial condition. Then

$${\hat{y}}^{\prime }(x)=2ax+3b{x}^2,\mathrm{ℰ}(x)=2ax+3b{x}^2-x+a{x}^2+b{x}^3.$$

Each method below gives two equations for the two unknowns $a$ and $b$.

1. The collocation method

Take $x_1=0.5$ and $x_2=1$. Then

$$\begin{array}{rl}\mathrm{ℰ}(0.5)& ={\hat{y}}^{\prime }(0.5)-0.5+\hat{y}(0.5)=1.25a+0.875b-0.5=0,\\ \mathrm{ℰ}(1)& ={\hat{y}}^{\prime }(1)-1+\hat{y}(1)=3a+4b-1=0.\end{array}$$

Solving these two equations gives

$$\hat{y}(x)=0.474x^2-0.105x^3,$$

an approximate solution in the sense of collocation.

If we know the exact solution, we can check the difference. Here $y(x)=x-1+e^{-x}$, and

$$|y(1/3)-\hat{y}(1/3)|\approx {10}^{-3},$$

so at $x=1/3$, for example, $\hat{y}$ is an approximate solution with accuracy about ${10}^{-3}$.

2. The method of moments

With two unknowns we use $1$ and $x$:

$${\int }_0^1(2ax+3b{x}^2-x+a{x}^2+b{x}^3)dx=0,{\int }_0^1(2ax+3b{x}^2-x+a{x}^2+b{x}^3)xdx=0.$$

These give

$$1.333a+1.25b-0.5=0,0.9167a+0.95b-0.333=0,$$

and so

$$\hat{y}(x)=0.483x^2-0.115x^3,$$

an approximate solution in the sense of the method of moments.

3. The Galerkin method

Solving the resulting system for $a,b$ gives

$$\hat{y}(x)=0.475x^2-0.107x^3,$$

the approximate solution in the sense of Galerkin.

4. The subdomain method

With the subintervals $[0,1/2]$ and $[1/2,1]$:

$${\int }_0^{1/2}\mathrm{ℰ}(x)dx=0,{\int }_{1/2}^1\mathrm{ℰ}(x)dx=0,$$

which gives

$$\hat{y}(x)=0.485x^2-0.118x^3$$

in the sense of the subdomain method.

5. The method of least squares

This gives

$$\hat{y}(x)=0.477x^2-0.110x^3$$

in the sense of least squares.

A boundary value problem

---

Galerkin Approximations in Two Dimensions

---

Natural Boundary Conditions

Integration identities

To handle natural boundary conditions we move derivatives from one function to another with the following identities.

1-D: integration by parts. On $\mathrm{\Omega }=[a,b]$ with boundary $\mathrm{\Gamma }=\{a,b\}$,

$${\int }_{\mathrm{\Omega }}udv=uv{|}_{\mathrm{\Gamma }}-{\int }_{\mathrm{\Omega }}vdu,$$

where $u(x),v(x)$ are integrable functions of $x$.

2-D: Green's theorem. For a region $\mathrm{\Omega }$ with closed boundary $\mathrm{\Gamma }$,

$${\oint }_{\mathrm{\Gamma }}\alpha dx+\beta dy={\iint }_{\mathrm{\Omega }}(\frac{\partial \beta }{\partial x}-\frac{\partial \alpha }{\partial y})d\mathrm{\Omega },$$

where $\alpha ,\beta$ are suitably defined differentiable and integrable functions.

3-D: the divergence theorem.

$${\int }_{\mathrm{\Gamma }}\overrightarrow{A}\cdot \overrightarrow{n}d\mathrm{\Gamma }={\int }_{\mathrm{\Omega }}\mathrm{\nabla }\cdot \overrightarrow{A}d\mathrm{\Omega },$$

where

$$\overrightarrow{A}=A_x\hat{i}+A_y\hat{j}+A_z\hat{k},\mathrm{\nabla }=\hat{i}\frac{\partial }{\partial x}+\hat{j}\frac{\partial }{\partial y}+\hat{k}\frac{\partial }{\partial z},\overrightarrow{n}=n_x\hat{i}+n_y\hat{j}+n_z\hat{k},\mathrm{\nabla }\cdot \overrightarrow{A}=\frac{\partial A_x}{\partial x}+\frac{\partial A_y}{\partial y}+\frac{\partial A_z}{\partial z}.$$

A natural boundary condition in one dimension

---

The General Weighted Residual Statement

Suppose

$$L\varphi +f=0\text{ in }\mathrm{\Omega },M\varphi +r=0\text{ on the boundary }\partial \mathrm{\Omega }=\mathrm{\Gamma },$$

where $L$ is a differential operator and $M$ is a boundary operator. For example,

$$L={\mathrm{\nabla }}^2=\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2},L={\mathrm{\nabla }}^4={\mathrm{\nabla }}^2({\mathrm{\nabla }}^2)=\frac{{\partial }^4}{\partial x^4}+2\frac{{\partial }^4}{\partial x^2\partial y^2}+\frac{{\partial }^4}{\partial y^4}.$$

Let $\hat{\varphi }=\psi +\sum _{m=1}^M{a}_m{N}_m$. The residuals are

$${\mathrm{ℰ}}_{\mathrm{\Omega }}=L\hat{\varphi }+f,{\mathrm{ℰ}}_{\mathrm{\Gamma }}=M\hat{\varphi }+r.$$

The weighted residual statement in general takes one of two forms:

$$\text{(1)}{\int }_{\mathrm{\Omega }}{\mathrm{ℰ}}_{\mathrm{\Omega }}wd\mathrm{\Omega }={\int }_{\mathrm{\Omega }}(L\hat{\varphi }+f)wd\mathrm{\Omega }=0,\text{(2)}{\int }_{\mathrm{\Omega }}{\mathrm{ℰ}}_{\mathrm{\Omega }}wd\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}{\mathrm{ℰ}}_{\mathrm{\Gamma }}\overline{w}d\mathrm{\Gamma }=0,$$

where $w$ and $\overline{w}$ are weight functions.

Form (1) with Galerkin weights $w_{\ell }=N_{\ell }$ gives

$$\sum _{m=1}^M{a}_m{\int }_{\mathrm{\Omega }}(L{N}_m)N_{\ell }d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}(L\psi )N_{\ell }d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}f{N}_{\ell }d\mathrm{\Omega }=0,\ell =1,\dots ,M,$$

that is, $[K]\{a\}=\{P\}$ with

$$K_{\ell m}={\int }_{\mathrm{\Omega }}{N}_{\ell }(L{N}_m)d\mathrm{\Omega },P_{\ell }=-{\int }_{\mathrm{\Omega }}{N}_{\ell }fd\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}{N}_{\ell }(L\psi )d\mathrm{\Omega }.$$

Form (2) with $w_{\ell }=N_{\ell }$ and ${\overline{w}}_{\ell }=-N_{\ell }$ (taking $\psi =0$) gives

$${\int }_{\mathrm{\Omega }}(L\hat{\varphi }+f)N_{\ell }d\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}(M\hat{\varphi }+r)N_{\ell }d\mathrm{\Gamma }=0,$$

$$\sum _{m=1}^M{a}_m{\int }_{\mathrm{\Omega }}(L{N}_m)N_{\ell }d\mathrm{\Omega }+{\int }_{\mathrm{\Omega }}f{N}_{\ell }d\mathrm{\Omega }-\sum _{m=1}^M{a}_m{\int }_{\mathrm{\Gamma }}(M{N}_m)N_{\ell }d\mathrm{\Gamma }-{\int }_{\mathrm{\Gamma }}r{N}_{\ell }d\mathrm{\Gamma }=0.$$

The terms multiplying $a_m$ form the matrix and the terms with $f$ and $r$ form the vector: $[K^{\prime }]\{a\}=\{P^{\prime }\}$ with

$$K_{\ell m}^{\prime }={\int }_{\mathrm{\Omega }}(L{N}_m)N_{\ell }d\mathrm{\Omega }-{\int }_{\mathrm{\Gamma }}(M{N}_m)N_{\ell }d\mathrm{\Gamma },P_{\ell }^{\prime }=-{\int }_{\mathrm{\Omega }}f{N}_{\ell }d\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}r{N}_{\ell }d\mathrm{\Gamma }.$$

---

The BEM for the Laplace and Poisson Equations

Consider

$$\{\begin{array}{ll}{\mathrm{\nabla }}^2\varphi =f& \text{in a region }\mathrm{\Omega }\subset {\mathrm{ℝ}}^2,f=f(x,y)\text{ a known function},\\ \varphi =\overline{\varphi }& \text{given on }{\mathrm{\Gamma }}_1\text{(forced, Dirichlet)},\\ q={\displaystyle \frac{\partial \varphi }{\partial n}}=\overline{q}& \text{given on }{\mathrm{\Gamma }}_2\text{(natural, Neumann)}.\end{array}$$

Let $\varphi$ be the exact solution (assume it is not available) and $\hat{\varphi }$ an approximate solution, with $\hat{q}=\partial \hat{\varphi }/\partial n$. Then the residuals are

$${\mathrm{ℰ}}_{\mathrm{\Omega }}={\mathrm{\nabla }}^2\hat{\varphi }-f,{\mathrm{ℰ}}_{{\mathrm{\Gamma }}_1}=\hat{\varphi }-\overline{\varphi },{\mathrm{ℰ}}_{{\mathrm{\Gamma }}_2}=\hat{q}-\overline{q}.$$

We use the Galerkin method to make the residuals as small as possible. The weighted residual statement is

$${\int }_{\mathrm{\Omega }}{\mathrm{ℰ}}_{\mathrm{\Omega }}{w}_{\ell }d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_1}{\mathrm{ℰ}}_{{\mathrm{\Gamma }}_1}{\overline{w}}_{\ell }d\mathrm{\Gamma }+{\int }_{{\mathrm{\Gamma }}_2}{\mathrm{ℰ}}_{{\mathrm{\Gamma }}_2}{\overline{\overline{w}}}_{\ell }d\mathrm{\Gamma }=0,$$

with the boundary weights

$${\overline{w}}_{\ell }=\frac{\partial w_{\ell }}{\partial n}\text{ on }{\mathrm{\Gamma }}_1,{\overline{\overline{w}}}_{\ell }=-w_{\ell }\text{ on }{\mathrm{\Gamma }}_2.$$

Substituting the residuals:

$${\int }_{\mathrm{\Omega }}({\mathrm{\nabla }}^2\hat{\varphi }-f)w_{\ell }d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_1}(\hat{\varphi }-\overline{\varphi })\frac{\partial w_{\ell }}{\partial n}d\mathrm{\Gamma }-{\int }_{{\mathrm{\Gamma }}_2}(\hat{q}-\overline{q})w_{\ell }d\mathrm{\Gamma }=0.$$

With $u=w_{\ell }$, $v=\hat{\varphi }$:

$${\int }_{{\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2}\frac{\partial \hat{\varphi }}{\partial n}{w}_{\ell }d\mathrm{\Gamma }-{\int }_{\mathrm{\Omega }}\mathrm{\nabla }\hat{\varphi }\cdot \mathrm{\nabla }{w}_{\ell }d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}f{w}_{\ell }d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_1}(\hat{\varphi }-\overline{\varphi })\frac{\partial w_{\ell }}{\partial n}d\mathrm{\Gamma }-{\int }_{{\mathrm{\Gamma }}_2}(\hat{q}-\overline{q})w_{\ell }d\mathrm{\Gamma }=0.$$

Apply Green's identity again, now with $u=\hat{\varphi }$, $v=w_{\ell }$:

$${\int }_{{\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2}\hat{q}{w}_{\ell }d\mathrm{\Gamma }-{\int }_{{\mathrm{\Gamma }}_1+{\mathrm{\Gamma }}_2}\hat{\varphi }\frac{\partial w_{\ell }}{\partial n}d\mathrm{\Gamma }+{\int }_{\mathrm{\Omega }}\hat{\varphi }{\mathrm{\nabla }}^2{w}_{\ell }d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}f{w}_{\ell }d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_1}(\hat{\varphi }-\overline{\varphi })\frac{\partial w_{\ell }}{\partial n}d\mathrm{\Gamma }-{\int }_{{\mathrm{\Gamma }}_2}(\hat{q}-\overline{q})w_{\ell }d\mathrm{\Gamma }=0.$$

On ${\mathrm{\Gamma }}_1$ the terms with $\hat{\varphi }$ cancel, and on ${\mathrm{\Gamma }}_2$ the terms with $\hat{q}$ cancel. Therefore

$$(\ast ){\int }_{\mathrm{\Omega }}\hat{\varphi }{\mathrm{\nabla }}^2{w}_{\ell }d\mathrm{\Omega }-{\int }_{\mathrm{\Omega }}f{w}_{\ell }d\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_1}(\hat{q}{w}_{\ell }-\overline{\varphi }\frac{\partial w_{\ell }}{\partial n})d\mathrm{\Gamma }+{\int }_{{\mathrm{\Gamma }}_2}(\overline{q}{w}_{\ell }-\hat{\varphi }\frac{\partial w_{\ell }}{\partial n})d\mathrm{\Gamma }=0.$$

The last two integrals are on the boundary; the first two are on the domain, and we need to convert them or get rid of them.

• Take $f=0$ to get rid of ${\int }_{\mathrm{\Omega }}f{w}_{\ell }d\mathrm{\Omega }$. (Later we will use the DRBEM to handle this inhomogeneous term with a domain integral.)
• We still have the domain integral ${\displaystyle {\int }_{\mathrm{\Omega }}\hat{\varphi }{\mathrm{\nabla }}^2{w}_{\ell }d\mathrm{\Omega }}$. To remove it, we choose $w_{\ell }$ so that ${\mathrm{\nabla }}^2{w}_{\ell }$ acts like the Dirac delta function.

The Dirac delta function

Thus in the weighted residual statement the multiplier ${\mathrm{\nabla }}^2{w}_{\ell }$ in the domain integral ${\int }_{\mathrm{\Omega }}\hat{\varphi }{\mathrm{\nabla }}^2{w}_{\ell }d\mathrm{\Omega }$ should be like the Dirac delta function:

$${\mathrm{\nabla }}^2{w}_{\ell }=\delta (|\overrightarrow{r}-{\overrightarrow{r}}_{\ell }|),$$

where ${\overrightarrow{r}}_{\ell }=(x_{\ell },y_{\ell })$ is the source point (like $x^{\prime }$ in $\delta (x,x^{\prime })$) and $\overrightarrow{r}=(x,y)$ is a general point, so $\overrightarrow{r}-{\overrightarrow{r}}_{\ell }$ is the position vector from the source point.

Back to the weighted residual statement $(\ast )$ with $f=0$:

$$c_{\ell }\hat{\varphi }({\overrightarrow{r}}_{\ell })+{\int }_{{\mathrm{\Gamma }}_2}(\overline{q}{w}_{\ell }-\hat{\varphi }\frac{\partial w_{\ell }}{\partial n})d\mathrm{\Gamma }+{\int }_{{\mathrm{\Gamma }}_1}(\hat{q}{w}_{\ell }-\overline{\varphi }\frac{\partial w_{\ell }}{\partial n})d\mathrm{\Gamma }=0,$$

where $c_{\ell }$ is a constant: $c_{\ell }=1$ for a source point inside $\mathrm{\Omega }$, and $c_{\ell }=\frac{1}{2}$ for a source point on a smooth part of the boundary (in one dimension, at an endpoint).

---

Fundamental Solutions

For the Laplace equation the fundamental solutions are:

$$\begin{array}{rl}\text{1-D:}& w_{\ell }=\frac{|x-x_{\ell }|}{2}=\frac{1}{2}[(x-x_{\ell })(H(x-x_{\ell })-H(x_{\ell }-x))],\\ \text{2-D:}& w_{\ell }=\frac{1}{2\pi }\ln |\overrightarrow{r}-{\overrightarrow{r}}_{\ell }|,\\ \text{3-D:}& w_{\ell }=\frac{-1}{4\pi |\overrightarrow{r}-{\overrightarrow{r}}_{\ell }|},\end{array}$$

where $H$ is the Heaviside function.

Deriving the fundamental solutions

Away from the source point, $w_{\ell }$ satisfies the Laplace equation, and by symmetry it depends only on the distance $r=|\overrightarrow{r}-{\overrightarrow{r}}_{\ell }|$.

2-D. The Laplacian in cylindrical coordinates is

$${\mathrm{\nabla }}^{2}u=\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^{2}u}{\partial {\theta }^2}+\frac{{\partial }^{2}u}{\partial z^2}=0.$$

Assume the axisymmetric case, so only $r$ varies:

$$\frac{d}{dr}\left(r\frac{du}{dr}\right)=0\Rightarrow r\frac{du}{dr}=A\Rightarrow \frac{du}{dr}=\frac{A}{r}\Rightarrow u=A\ln r+B.$$

3-D. Use spherical coordinates, with $\varphi$ the angle from the $z$-axis and $\theta$ the angle in the $xy$-plane:

$${\mathrm{\nabla }}^{2}u=\frac{1}{r^2}[\frac{\partial }{\partial r}\left(r^2\frac{\partial u}{\partial r}\right)+\frac{1}{{\sin }^2\varphi }\frac{{\partial }^{2}u}{\partial {\theta }^2}+\frac{1}{\sin \varphi }\frac{\partial }{\partial \varphi }(\sin \varphi \frac{\partial u}{\partial \varphi })]=0.$$

Assuming axisymmetry we have only the variable $r$:

$$\frac{d}{dr}\left(r^2\frac{du}{dr}\right)=0\Rightarrow r^2\frac{du}{dr}=A\Rightarrow \frac{du}{dr}=\frac{A}{r^2}\Rightarrow u=C\frac{1}{r}+B.$$

Computing the constants. Take a simple domain $D_{\epsilon }$ around the point $x^{\prime }={\overrightarrow{r}}_{\ell }=(x_{\ell },y_{\ell })$: an interval of length $2\epsilon$ in 1-D, a disc of radius $\epsilon$ in 2-D, a ball of radius $\epsilon$ in 3-D. Integrating ${\mathrm{\nabla }}^2{w}_{\ell }=\delta$ over $D_{\epsilon }$,

$${\int }_{D_{\epsilon }}{\mathrm{\nabla }}^2{w}_{\ell }dV={\int }_{D_{\epsilon }}\delta (x,x^{\prime })dV=1.$$

Apply the divergence theorem:

$${\int }_{\partial D_{\epsilon }}\frac{\partial w_{\ell }}{\partial n}d\mathrm{\Gamma }=1.$$

On $\partial D_{\epsilon }$ the normal derivative is $\partial w_{\ell }/\partial r$ evaluated at $r=\epsilon$, which is constant, so the integral is this value times the size of $\partial D_{\epsilon }$.

• In 3-D: ${\displaystyle \frac{d{w}_{\ell }}{dr}}{|}_{r=\epsilon }=-{\displaystyle \frac{C}{{\epsilon }^2}}$ and the sphere has area $4\pi {\epsilon }^2$, so $-{\displaystyle \frac{C}{{\epsilon }^2}}\cdot 4\pi {\epsilon }^2=1$, i.e. $C=-{\displaystyle \frac{1}{4\pi }}$.
• In 2-D: ${\displaystyle \frac{d{w}_{\ell }}{dr}}{|}_{r=\epsilon }={\displaystyle \frac{A}{\epsilon }}$ and the circle has length $2\pi \epsilon$, so ${\displaystyle \frac{A}{\epsilon }}\cdot 2\pi \epsilon =1$, i.e. $A={\displaystyle \frac{1}{2\pi }}$.
• In 1-D: $w_{\ell }=A|x-x_{\ell }|$, and the two endpoints contribute $A+A=1$, so $A={\displaystyle \frac{1}{2}}$.

In both 2-D and 3-D the fundamental solutions are, up to a constant,

$$w_{\ell }=\frac{1}{2\pi }\ln r+B,w_{\ell }=-\frac{1}{4\pi r}+B,$$

and $B=0$ is taken from the boundary condition on $w_{\ell }$.

---

A One-Dimensional Example

---

Flashcards

One card for every definition in these notes. Tap or click a card (or press the space bar) to see the definition; use the arrows to move between cards.