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$: egin{equation} label{eq:harmonic-function} riangle u = 0 quad (x in mathbb{R}^d). end{equation} Let $gamma(x, y)$ be the fundamental solution for the free space such that egin{equation} label{eq:laplace-equation} - riangle_x gamma(x, y) = delta(x - y) quad (x, y in mathbb{R}^d). end{equation} It has the following formulation: egin{equation} label{eq:fundamental-solution} gamma(x, y) = egin{cases} -frac{1}{2pi}lnlvert x - y vert & (d = 2) \ frac{lvert x - y vert^{2-d}}{(d-2)omega_d} & (d > 2) end{cases}, end{equation} where $omega_d = frac{2pi^{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 = pdiffOmega$. Let $vect{F} = psi ablavarphi - varphi ablapsi$ and apply the Gauss divergence theorem, we have the famous Green's 2nd identity as below: egin{equation} label{eq:green-2nd-identity} int_{Omega} left( psi rianglevarphi - varphi rianglepsi ight) intd V = int_{Gamma} left( psi frac{pdiffvarphi}{pdiff ormvect} - varphi frac{pdiff psi}{pdiff ormvect} ight) intd S, end{equation} where $ ormvect$ 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 egin{equation} label{eq:green-2nd-identity-with-fundamental-solution} int_{Omega} left( gamma(x,y) riangle_x u(x) - u(x) riangle_xgamma(x,y) ight) intd V(x) = int_{Gamma} left( gamma(x,y) frac{pdiff u(x)}{pdiff ormvect(x)} - u(x) frac{pdiff gamma(x,y)}{pdiff ormvect(x)} ight) 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 ormvect(x)} intd S(x) - int_{Gamma} u(x) frac{pdiff gamma(x,y)}{pdiff ormvect(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 egin{align} label{eq:fundamental-solution-symmetry} gamma(x,y) &= gamma(y,x) \ frac{pdiffgamma(y,x)}{pdiff ormvect(y)} = K^{*}(y,x) &= K(x,y) = frac{pdiffgamma(x,y)}{pdiff ormvect(y)}, end{align} after swappping the variables $x$ and $y$, we have the representation formula for the interior $Int(Omega)$ of $Omega$ as below: egin{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 ormvect(y)}$ and $K(x,y) = frac{pdiffgamma(x,y)}{pdiff ormvect(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 ackslash overline{Omega}$ of $Omega$, a representation formula with the same formulation can be obtained as long as we assume that when $abs{x} ightarrow 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 egin{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 ormvect'(y)}$ and $K'(x,y) = frac{pdiffgamma(x,y)}{pdiff ormvect'(y)}$, where $ ormvect'$ is the unit outward normal vector with respect to the domain $Omega'$, which has opposite direction compared to $ ormvect$.
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 ackslash Gamma) $$ with a charge density $phi$. This equation has the following two cases:
- 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 $ ormvect$ be the unit outward normal vector of $Omega$, the function $$ u_t(x) = u(x + t ormvect(x)) quad (x in Gamma) $$ converges uniformly to $u_-$ when $t ightarrow 0^-$ and to $u_+$ when $t ightarrow 0^+$, where egin{equation} egin{aligned} u_{-}(x) &= -frac{1}{2} phi(x) + T_Kphi = -frac{1}{2} phi(x) + int_{Gamma} K(x, y) phi(y) intd o(y) \ u_{+}(x) &= frac{1}{2} phi(x) + T_Kphi = 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} egin{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 = egin{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 egin{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' = egin{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 $ ormvect$ with respect to $Omega$ in eqref{eq:interior-representation-formula-behavior}, we have egin{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$ egin{equation} label{eq:double-layer-potential-jump} int_{Gamma} K(x,y) u(y)igvert_{Omega'} intd S(y) - int_{Gamma} K(x,y) u(y)igvert_{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.