This article introduces how to derive the representation formula used in BEM from Green‘s identity.
Interior and exterior representation formulas
Let $u$ be a harmonic function in the free space $\mathbb{R}^d$: \begin{equation} \label{eq:harmonic-function} \triangle u = 0 \quad (x \in \mathbb{R}^d). \end{equation} Let $\gamma(x, y)$ be the fundamental solution for the free space such that \begin{equation} \label{eq:laplace-equation} -\triangle_x \gamma(x, y) = \delta(x - y) \quad (x, y \in \mathbb{R}^d). \end{equation} It has the following formulation: \begin{equation} \label{eq:fundamental-solution} \gamma(x, y) = \begin{cases} -\frac{1}{2\pi}\ln\lvert x - y \rvert & (d = 2) \\ \frac{\lvert x - y \rvert^{2-d}}{(d-2)\omega_d} & (d > 2) \end{cases}, \end{equation} where $\omega_d = \frac{2\pi^{d/2}}{\Gamma(d/2)}$, $x$ is the field point and $y$ is the source point. Let $\psi$ and $\varphi$ be two functions having 2nd order derivatives in a bounded domain $\Omega$ in $\mathbb{R}^d$ with its boundary $\Gamma = \pdiff\Omega$. Let $\vect{F} = \psi\nabla\varphi - \varphi\nabla\psi$ and apply the Gauss divergence theorem, we have the famous Green‘s 2nd identity as below: \begin{equation} \label{eq:green-2nd-identity} \int_{\Omega} \left( \psi\triangle\varphi - \varphi\triangle\psi \right) \intd V = \int_{\Gamma} \left( \psi \frac{\pdiff\varphi}{\pdiff \normvect} - \varphi \frac{\pdiff \psi}{\pdiff \normvect} \right) \intd S, \end{equation} where $\normvect$ is the unit outward normal vector with respect to domain $\Omega$, which points from interior to exterior. By replacing $\psi$ with $\gamma(x,y)$ and $\varphi$ with $u(x)$, and performing integration and differentiation with respect to the variable $x$, we have \begin{equation} \label{eq:green-2nd-identity-with-fundamental-solution} \int_{\Omega} \left( \gamma(x,y)\triangle_x u(x) - u(x)\triangle_x\gamma(x,y) \right) \intd V(x) = \int_{\Gamma} \left( \gamma(x,y) \frac{\pdiff u(x)}{\pdiff \normvect(x)} - u(x) \frac{\pdiff \gamma(x,y)}{\pdiff \normvect(x)} \right) \intd S(x). \end{equation} After substituting \eqref{eq:harmonic-function} and \eqref{eq:laplace-equation}, we have $$ u(y) = \int_{\Gamma} \gamma(x,y) \frac{\pdiff u(x)}{\pdiff \normvect(x)} \intd S(x) - \int_{\Gamma} u(x) \frac{\pdiff \gamma(x,y)}{\pdiff \normvect(x)} \intd S(x) \quad (y \in \Int(\Omega)). $$ where $\Int(\Omega)$ is the interior of $\Omega$. Due to the symmetric property of the fundamental solution \begin{align} \label{eq:fundamental-solution-symmetry} \gamma(x,y) &= \gamma(y,x) \\ \frac{\pdiff\gamma(y,x)}{\pdiff \normvect(y)} = K^{*}(y,x) &= K(x,y) = \frac{\pdiff\gamma(x,y)}{\pdiff\normvect(y)}, \end{align} after swappping the variables $x$ and $y$, we have the representation formula for the interior $\Int(\Omega)$ of $\Omega$ as below: \begin{equation} \label{eq:interior-representation-formula} u(x) = \int_{\Gamma} \gamma(x,y) \psi(y) \intd S(y) - \int_{\Gamma} K(x,y) u(y) \intd S(y) \quad (x \in \Int(\Omega)), \end{equation} where $\psi(y) = \frac{\pdiff u(y)}{\pdiff \normvect(y)}$ and $K(x,y) = \frac{\pdiff\gamma(x,y)}{\pdiff\normvect(y)}$. The first term in the above equation is the single layer potential, while the second term is the double layer potential.
Remark It can be seen that the interior representation formula in equation \eqref{eq:interior-representation-formula} has the same formulation as that derived from the direct method.
For the exterior $\Omega‘ = \mathbb{R}^d \backslash \overline{\Omega}$ of $\Omega$, a representation formula with the same formulation can be obtained as long as we assume that when $\abs{x} \rightarrow \infty$, both $\gamma(x,y)$ and $K(x,y)$ approach to zero, so that the integration on infinite boundary has no contribution. Therefore, the representation formula for the exterior of $\Omega$ is \begin{equation} \label{eq:exterior-representation-formula} u(x) = \int_{\Gamma} \gamma(x,y) \psi‘(y) \intd S(y) - \int_{\Gamma} K‘(x,y) u(y) \intd S(y) \quad (x \in \Int(\Omega‘)). \end{equation} Here $\psi‘(y) = \frac{\pdiff u(y)}{\pdiff \normvect‘(y)}$ and $K‘(x,y) = \frac{\pdiff\gamma(x,y)}{\pdiff\normvect‘(y)}$, where $\normvect‘$ is the unit outward normal vector with respect to the domain $\Omega‘$, which has opposite direction compared to $\normvect$.
Representation formula at the boundary $\Gamma$
It is well known that the single layer potential in equation \eqref{eq:interior-representation-formula} or \eqref{eq:exterior-representation-formula} is continuous across the boundary $\Gamma$, while the double layer potential has a jump, which is governed by the following theorem.
Theorem (Boundary limit of double layer potential) Let $\phi \in C(\Gamma)$ and $u$ be the double layer potential $$ u(x) = \int_{\Gamma} K(x,y) \phi(y) \intd S(y) \quad (x \in \mathbb{R}^d \backslash \Gamma) $$ with a charge density $\phi$. This equation has the following two cases: \alert{This needs to be further clarified!}
- interior representation formula: $$ u(x) = \int_{\Gamma} K(x,y) \phi(y) \intd S(y) \quad (x \in \Omega) $$
- exterior representation formula: $$ u(x) = \int_{\Gamma} K‘(x,y) \phi(y) \intd S(y) \quad (x \in \Omega‘) $$
Then the restrictions of $u$ to $\Omega$ and $\Omega‘$ both have continuous extensions to $\overline{\Omega}$ and $\overline{\Omega‘}$ respectively. Let $t \in \mathbb{R}$ and $\normvect$ be the unit outward normal vector of $\Omega$, the function $$ u_t(x) = u(x + t \normvect(x)) \quad (x \in \Gamma) $$ converges uniformly to $u_-$ when $t \rightarrow 0^-$ and to $u_+$ when $t \rightarrow 0^+$, where \begin{equation} \begin{aligned} u_{-}(x) &= -\frac{1}{2} \phi(x) + T_K\phi = -\frac{1}{2} \phi(x) + \int_{\Gamma} K(x, y) \phi(y) \intd o(y) \\ u_{+}(x) &= \frac{1}{2} \phi(x) + T_K\phi = \frac{1}{2} \phi(x) + \int_{\Gamma} K(x, y) \phi(y) \intd o(y) \end{aligned} \quad (x \in \Gamma). \end{equation}
Representation formula outside the domain $\Omega$
For the interior representation formula \eqref{eq:interior-representation-formula}, when the variable $x$ is outside the domain $\Omega$, $u$ evaluates to zero. This is because according to equation \eqref{eq:green-2nd-identity-with-fundamental-solution}, before swapping $x$ and $y$, when the variable $y$ is outside $\Omega$, the Dirac function $\Delta_x \gamma(x,y) = -\delta(x - y)$ evaluates to zero. Similarly, for the exterior representation formula \eqref{eq:exterior-representation-formula}, when the variable $x$ is outside the domain $\Omega‘$, $u$ also evaluates to zero.
Summary of representation formulas‘ behavior in $\mathbb{R}^d$
By summarizing previous results, we can conclude that for the interior representation formula \eqref{eq:interior-representation-formula} \begin{equation} \label{eq:interior-representation-formula-behavior} \int_{\Gamma} \gamma(x,y) \psi(y) \intd S(y) - \int_{\Gamma} K(x,y) u(y) \intd S(y) = cu(x) \end{equation} where $$ c = \begin{cases} 1 & x \in \Int(\Omega) \\ \frac{1}{2} & x \in \Gamma \\ 0 & x \in \Int(\Omega‘) \end{cases} $$ Similarly for the exterior representation formula \eqref{eq:exterior-representation-formula} we have \begin{equation} \label{eq:exterior-representation-formula-behavior} \int_{\Gamma} \gamma(x,y) \psi‘(y) \intd S(y) - \int_{\Gamma} K‘(x,y) u(y) \intd S(y) = c‘u(x) \end{equation} where $$ c‘ = \begin{cases} 1 & x \in \Int(\Omega‘) \\ \frac{1}{2} & x \in \Gamma \\ 0 & x \in \Int(\Omega) \end{cases} $$ If we also use the normal vector $\normvect$ with respect to $\Omega$ in \eqref{eq:interior-representation-formula-behavior}, we have \begin{equation} \label{eq:interior-representation-formula-with-normvect} -\int_{\Gamma} \gamma(x,y) \psi(y) \intd S(y) + \int_{\Gamma} K(x,y) u(y) \intd S(y) = c‘u(x). \end{equation} It should be noted that although the left hand sides of \eqref{eq:interior-representation-formula-behavior} and \eqref{eq:interior-representation-formula-with-normvect} have the same form with opposite signs, they do not cancel with other because the limiting values of the double layer charge density $u$ used in the integral are approached to $\Gamma$ from interior and exterior respectively. Therefore, although the single layer potential is continuous across the boundary $\Gamma$, the double layer potential has a jump. Then we have the following jump behavior for the double layer potential at $\Gamma$ \begin{equation} \label{eq:double-layer-potential-jump} \int_{\Gamma} K(x,y) u(y)\big\vert_{\Omega‘} \intd S(y) - \int_{\Gamma} K(x,y) u(y)\big\vert_{\Omega} \intd S(y) = u(x) \quad (x \in \Gamma), \end{equation} which is consistent with $u_+ - u_- = \phi$ derived from the theorem for the boundary limit of double layer potential.