% % 2009-04-25
\documentclass[12pt]{amsart}
%\usepackage{a4wide}
%\addtolength{\hoffset}{1.45cm}
%\addtolength{\textwidth}{-2.9cm}
%\addtolength{\voffset}{1.2cm}
%\addtolen\texttt{}gth{\textheight}{-2.4cm}
%\usepackage{breqn}
%\usepackage{showkeys}

\usepackage{amsthm,amssymb,amsmath}
\usepackage{geometry}
\usepackage{graphicx}
\usepackage{enumerate}
\usepackage{float}

%\usepackage[active]{srcltx}
\numberwithin{equation}{section}
%\usepackage{pstricks,pst-plot}
\newtheorem{example}{example}[section]

\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\Span}{span}
\DeclareMathOperator{\diam}{diam}
\DeclareMathOperator*{\essup}{ess-sup}
\DeclareMathOperator{\Tr}{Tr}

\newtheorem{thm}{Theorem}
\numberwithin{thm}{section}

\newtheorem{lem}[thm]{Lemma}
\newtheorem{coro}[thm]{Corollary}
\newtheorem{propo}[thm]{Proposition}
\theoremstyle{definition}
%\newtheorem{defi}[thm]{Definition}
\newtheorem{ex}[thm]{Example}
\newtheorem{rem}[thm]{Remark}

\newcommand{\inner}[3][]{\langle #2,#3\rangle_{#1}}
\newcommand{\Inner}[3][]{\Big\langle#2,#3\Big\rangle_{#1}}
\newcommand{\norm}[2][]{\|{#2}\|_{{#1}}}
\newcommand{\snorm}[2][]{|{#2}|_{{#1}}}
\newcommand{\Snorm}[2][]{\Big|{#2}\Big|_{{#1}}}
\newcommand{\Norm}[2][]{\Big\|{#2}\Big\|_{{#1}}}
\newcommand{\abs}[1]{|#1|}
\newcommand{\Abs}[1]{\Big|#1\Big|}
\newcommand{\rD}{{\mathrm D}}
\newcommand{\TT}{{\mathcal T}}
\newcommand{\EE}{{\mathbf E}}
\newcommand{\dd}{ { \mathrm{d}} }
\newcommand{\ee}{ { \mathrm{e}} }
\newcommand{\ii}{ { \mathrm{i}} }

\newcommand{\pt}{\partial}

\newcommand{\cT}{{ \mathcal T}}
\newcommand{\cG}{{ \mathcal G}}
\newcommand{\cV}{{ \mathcal V}}
\newcommand{\cW}{{ \mathcal W}}
\newcommand{\bbT}{{ \mathbb T}}
\newcommand{\cZ}{{ \mathcal Z}}
\newcommand{\cR}{{ \mathcal R}}
\newcommand{\cK}{{ \mathcal K}}
\newcommand{\cA}{{ \mathcal A}}
\newcommand{\cP}{{ \mathcal P}}
\newcommand{\cC}{{ \mathcal C}}
\newcommand{\cH}{{ \mathcal H}}
\newcommand{\cE}{{ \mathcal E}}
\newcommand{\cB}{{ \mathcal B}}
\newcommand{\cL}{{ \mathcal L}}
\newcommand{\cF}{{ \mathcal F}}
\newcommand{\cD}{{ \mathcal D}}
\newcommand{\cN}{{ \mathcal N}}
\newcommand{\cO}{{ \mathcal O}}
\newcommand{\cM}{{ \mathcal M}}
\newcommand{\IR}{{ \mathbf R}}
\newcommand{\IC}{{ \mathbf C}}
\newcommand{\IT}{{ \mathbf T}}
\newcommand{\IQ}{{ \mathbf Q}}
\newcommand{\IP}{{ \mathbf P}}
\newcommand{\IN}{{ \mathbf N}}
\newcommand{\IE}{{ \mathbf E}}
\newcommand{\IInt}{{ \textrm Int}}
\newcommand{\HS}{\mathrm{HS}}


\title{Boundary element method for the Helholtz equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\author[A.~Mesforush]{Ali Mesforush}

\address{
  Department of Mathematical Sciences,
  Shahrood University,
  Shahrood, 
  Iran}

\email{ali.mesforush@shahroodut.ac.ir}
\urladdr{http://shahroodut.ac.ir/fa/as/?id=S138}

\author[Z.~Mehraban]{Zohreh Mehraban}

\address{
 Department of Mathematical Sciences,
  Shahrood University,
  Shahrood, 
  Iran}

\email{zohre.mehraban@ymail.com}
%\urladdr{http://www.math.chalmers.se/~mesforus}


\keywords{Laplace equation, Helmholtz equation, boundary integral equation, discretization, boundary element method}

\subjclass{35A08,65M38,65N38,74S15}

\begin{document}
\begin{abstract}
The boundary element methos is applied to the Helholtltz equation. To this end we discretized the Helmholtz equation over the boundary $\Omega$ and conclude a system of equations. By applying the boundary condition we get a new system which gives the solution of the equation.
\end{abstract}
\maketitle
%\textbf{Keywords:}
%Laplace equation, Helmholtz equation, boundary integral equation, discretization, boundary element method
\section{Introduction and preliminaries }
Numerical approximation for solving a differential equation has a wide range of application in engineering and mathematics. Some of these methods are: finite difference methods, finite element methods and boundary element method. In these methods we usually discritize the equation and apply the governing equations on the suitable mesh. The most important thing in the numerical approximation is determination of the size of mesh such that the approximation solution has the suitable accuracy.


The boundary element method is a powerful tool for numerical approximation for a differential equation. For the first time, in 1872, Betti presented this method in elasticity theory\cite{betti}. In 1903 Fredholm extended this method\cite{fred}, but in 1977 the boundary element method appears in some publications  of Banerjee, Butterfield, Brebbia and Dominguez.\cite{banerjee,doming}\\
In this paper we want to find an approximate solution for the Helmholtz equation by the boundary element method.  This  paper organized in three sections: in  the first section , we introduce some preliminaries and we describe the structure of the boundary element method. In the second section, we applied boundary element method for Laplace equation which is a special case of helmholtz equation. Finally, in the third section we  approximat the boundary element solution of Helmholtz equation.


%\section{preliminaries}
In boundary element method we use three important things.\cite{boundary}
\begin{itemize}
\item[1-]
{\it\textit{The Gauss divergence theorem:}}
If $\Omega \subset \mathbb{R}^d$ be a bounded domain and $w=(w_1,\cdots,w_d)\in C^1(\overline{\Omega})^d$ be a vector-valued function, then for $x\in \mathbb{R}^d$ we have
\begin{equation*}
\int_{\Omega}\nabla\cdot w\,\dd v=\int_{\Gamma}w\cdot n\,\dd s,
\end{equation*}
where $n=(n_1,\cdots,n_d)$ is the outward unit normal to boundary $\Gamma$.
\item[2-]
{\it\textit{The second Green identity:}}
Let $\Omega \subset \mathbb{R}^d$ be a bounded domain and $v, w\in C^2(\overline{\Omega})$, then
\begin{equation*}
\int_\Omega (w\nabla^2u -u\nabla^2 w)\dd v=\int_{\Gamma}(w\dfrac{\partial u}{\partial n}-u\dfrac{\partial w}{\partial n})\dd s.
\end{equation*}
\item[3-]
{\it\textit{The Dirac delta function:}}
The Dirac delta function defined as
\begin{equation*}
\delta (X-X_0)=
\begin{cases}
\infty,\qquad &X=X_0,\\
0, &X\neq X_0,
\end{cases}
\end{equation*} 
and has two following properties:
\begin{equation*}
\int_{-\infty}^{\infty}\delta (X-X_0)\,\dd X=1, \,\,
\int_{-\infty}^{\infty}f(X) \delta (X-X_0)\,\dd X=f(X_0).
\end{equation*}
\end{itemize}
%%%%%%%%%%%%%%%%%%%%%%

The structure of boundary element method can be decribed as follow.

 
%\section{Boundary element method Structure}
Consider  the following boundary element value problem 
\begin{align}\label{am1}
\begin{cases}
\mathcal{A}u:=Au_{xx}+2Bu_{xy}+Cu_{yy}+Du=f,& \mathrm{ in}\, \Omega,\\
u=\bar{u},& \mathrm{on}\, \Gamma_1,\\
q=\dfrac{\partial u}{\partial n}=\bar{q},& \mathrm{on}\,\Gamma_2,
\end{cases}
\end{align}
where the coefficients $A, B, C, D $ are constant and  $\Gamma_1$, $\Gamma_2$ are boundaries of domain $\Omega$, $q$ is the derivative of potential function in direction of $n$ and $\bar{u}$, $\bar{q}$
 are given values of the flux and potential function on the boundary. In boundary element method we multiply the first equation of \eqref{am1} by weight function $w\in C^2(\overline{\Omega})$ and integrate over $\Omega$ to get the integral form
 \begin{equation}\label{ta1}
 \int_{\Omega}(\mathcal{A}u)w\dd v=\int_\Omega fw\dd v.
 \end{equation}
 by choosing the fundamental solution $\mathcal{A}w=-\delta(X-X_0)$ as the weight function and using the second Green identity, we can convert \eqref{ta1} to a continuous integral equation over the boundary. For applying the boundary element method we  discretize the problem over the boundary of $\Omega$. To this end we consider $\partial\Omega=\cup_{j=1}^{N}\partial\Omega_j$ and conclude the system of equations $HU=GQ$. If we apply the boundary conditions, we get a new system of equations as $AX=B$, which gives the flux and potential $u$ on the boundary.\\
In this paper we want to apply the bounday element method for the Helmholtz equation. In \eqref{am1} if we consider $A=1, B=0, C=1, D=k^2$ and by considering $\nabla^2u=u_{xx}+u_{yy}$, we get the Helmholtz equation:
\begin{equation}\label{ma3}
\begin{cases}
\mathcal{A}u:=\nabla^2 u+k^2u=0, &\mathrm{ in}\,\Omega\subset\mathbb{R}^2,\\
u=\bar{u},& \mathrm{on}\,\Gamma_1,\\
q=\dfrac{\partial u}{\partial n}=\bar{q},&\mathrm{on}\, \Gamma_2,
\end{cases}
\end{equation}
where $k$ is the wave number.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Boundary element method for the laplace problem}
In this section we apply the boundary element method to the Laplace equation which is the special case of Helmholtz equation. 
In equation \eqref{ma3} if we set $k=0$, we get the laplace equation,
\begin{equation}\label{ma4}
\begin{cases}
\nabla^2 u=0, &\mathrm{in}\,\Omega\subset\mathbb{R}^2,\\
u=\bar{u},& \mathrm{on}\,\Gamma_1,\\
q=\dfrac{\partial u}{\partial n}=\bar{q},&\mathrm{on}\, \Gamma_2.
\end{cases}
\end{equation}
Consider weight function $w$ such that it satisfies 
\begin{equation}\label{zar}
\nabla^2 w+\delta(X-X_0)=0,
\end{equation}
 where $X_0=(\xi,\eta)$ is a singular point and $\delta$ is the Dirac delta function. Considering the polar form $\nabla^2 w =\dfrac{1}{r}\dfrac{\partial}{\partial r}(r\dfrac{\partial w}{\partial r})$ where $r=|X-X_0|=\sqrt{(x-\xi)^2+(y-\eta)^2}$. For $r>0, \delta(X-X_0)=0 $, thus
\begin{equation*}
\dfrac{1}{r}\dfrac{\partial}{\partial r}(r\dfrac{\partial w}{\partial r})=0,
\end{equation*}
where $A$ and $B$ are arbitrary constants. Since we look for a particular solution, we
may set $B=0$.
Which gives after integrating twice
\[w= A\ln r + B.\]
 By considering a neighborhood of $X_0$ and using the integral property  for dirac delta function, we obtain $A=\dfrac{1}{2\pi}$. Hence, the fundamental solution for the Laplac equation  becomes \cite{brebbia}
\[w=\dfrac{1}{2\pi}\ln \left(\dfrac{1}{r}\right).\]

Also we note that
\begin{equation*}
\dfrac{\partial w}{\partial n}=\dfrac{\partial w}{\partial r}\dfrac{\partial r}{\partial n}=-\dfrac{1}{2\pi r}({\nabla} r\cdot n)=-\dfrac{1}{2\pi r^2}(r_xn_x+r_yn_y),
\end{equation*}
where $r_x=x-\xi$ and $r_y=y-\eta$. 

 Let $X_0\in \Omega$ and consider a neighborhood of $X_0$ with radius $\varepsilon$, namely $\Omega_\varepsilon$ and consider the new domain $\Omega-\Omega_{\varepsilon}$ with boundary $\partial\Omega\cup\partial\Omega_\varepsilon$. In the new domain  by letting  $q=\dfrac{\partial u}{\partial n}$ and $q^*=\dfrac{\partial w}{\partial n}$ the Green formula can be written as
\begin{align*}
\lim_{\varepsilon\rightarrow 0}\int_{\Omega -\Omega_\varepsilon}(u\nabla^2 w-w\nabla^2 u)\dd v=\lim_{\varepsilon\rightarrow 0}&\int_{\partial\Omega }(uq^*-wq)\dd s\\
&+\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon }(uq^*-wq)\dd s,
\end{align*}
by using \eqref{zar} and the definition of the Dirac delta and $X_0\notin\Omega -\Omega_\varepsilon$ obtain
\[\lim_{\varepsilon\rightarrow 0}\int_{\Omega -\Omega_\varepsilon}u\nabla^2 w\dd v=\lim_{\varepsilon\rightarrow 0}\left(-\int_{\Omega -\Omega_\varepsilon}u\delta(X-X_0)\dd v\right)=0,
\]
by using the Laplac equation \eqref{ma4} obtain
\begin{equation*}
\lim_{\varepsilon\rightarrow 0}\int_{\Omega -\Omega_\varepsilon}(u\nabla^2 w-w\nabla^2u)\dd v=0.
\end{equation*}
In other hand,  \cite{wrobel}
\begin{align*}
%\begin{split}
\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon}(uq^*-wq)\dd s &=\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon}\left(u-u(X_0)+u(X_0)\right)q^*\dd s\\
& \quad -\lim_{\varepsilon\rightarrow 0}\left(-\dfrac{\ln \varepsilon}{2\pi}\int_{\partial\Omega_\varepsilon}\dfrac{\partial u}{\partial n}\dd s\right)\\
&=\lim_{\varepsilon\rightarrow 0}\left(u(X_0)\int_0^{2\pi}\dfrac{1}{2\pi}\dd\theta\right)+\lim_{\varepsilon\rightarrow 0}\left(-\dfrac{\varepsilon\ln\varepsilon}{2\pi}\dfrac{\partial u(X_0)}{\partial n}\int_0^{2\pi}\dd\theta\right)\\
&=u(X_0).
%\end{split}
\end{align*}
So the boundary integral equation in $\Omega$ is:
\begin{equation*}
u(X_0)=\int_{\partial\Omega }(wq-uq^*)\dd s,
\end{equation*}
which means that for given $u$, $q$ on $\partial\Omega$, the values of $u$ can be determinded for the interior point of the domain.\\
Let $X_0\in\partial\Omega$ and consider a neighborhood of $X_0$, namely $\Omega_\varepsilon^i$ and consider the new domain $\Omega-\Omega_{\varepsilon}^i$ with boundary $\partial\Omega_\varepsilon^i\cup(\partial\Omega-\partial\Omega_\varepsilon)$. In the new domain the Green formula can be written as,
\begin{align*}
\lim_{\varepsilon\rightarrow 0}\int_{\Omega -\Omega_\varepsilon^i}(u\nabla^2 w-w\nabla^2 u)\dd v & =\lim_{\varepsilon\rightarrow 0}\int_{ (\partial\Omega-\partial\Omega_\varepsilon)}(uq^*-wq)\dd s\\
&\quad +\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon^i}(uq^*-wq)\dd s.
\end{align*}
Since the governing equation is Laplace equation and the $X_0$ is out of  the domain $\Omega -\Omega_\varepsilon^i$,we have:
\begin{equation*}
\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega -\partial\Omega_\varepsilon}(u\nabla^2 w-w\nabla^2u)\dd v=0.
\end{equation*}
In other hand
\begin{equation*}
\begin{split}
\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon^i}(uq^*-wq)\dd s &=\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon}\left(u-u(X_0)+u(X_0)\right)q^*\dd s\\
& \quad -\lim_{\varepsilon\rightarrow 0}\left(-\dfrac{\ln \varepsilon}{2\pi}\int_{\partial\Omega_\varepsilon}\dfrac{\partial u}{\partial n}\dd s\right)\\
&=\lim_{\varepsilon\rightarrow 0}\left(u(X_0)\int_{\partial\Omega_\varepsilon^i}\dfrac{1}{2\pi}\dd\theta\right)+\lim_{\varepsilon\rightarrow 0}\left(-\dfrac{\varepsilon\ln\varepsilon}{2\pi}\dfrac{\partial u(X_0)}{\partial n}\int_0^{2\pi}\dd\theta\right)\\
&=C(X_0)u(X_0),
\end{split}
\end{equation*}
where $C(X_0)$ is a geometry-dependent parameter. It varies between zero and unity and equals $\frac{1}{2}$. If the point is on a smooth part of the boundary at edges, it is related to the angle of the joining surfaces and equals $C(X_0)=\dfrac{\text{internal angle}}{2\pi}$. Also
\begin{equation*}
\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega_\varepsilon}(uq^*-wq)\dd s =\lim_{\varepsilon\rightarrow 0}\int_{\partial\Omega}(uq^*-wq)\dd s.
\end{equation*}
So the boundary integral equation in $\partial\Omega$ is:
\begin{equation*}
C(\alpha)u(X_0)+\int_{\partial\Omega }uq^*\dd s=\int_{\partial\Omega } wq\dd s.
\end{equation*}
To apply the boundary element method, the boundary is divided into small
elements as $\partial\Omega=\cup_{j=1}^{N}\partial\Omega_j$ where $N$ represents the number of elements that form the boundary. The mesh here is generated only on the boundary and it is more flexible than in the case of finite element mesh. The integrals are then expressed for each element and summed up in order to evaluate the integral on the whole boundary. Considering constant element discretization and the source points situated on a smooth boundary, the boundary  integral equation can be written as
\begin{equation*}
\dfrac{1}{2}u(X_0)+\int_{\cup_{j=1}^{N}\partial\Omega_j}uq^*\dd s=\int_{\cup_{j=1}^{N}\partial\Omega_j } wq\dd s.
\end{equation*}
where by definition
\begin{equation*}
H_{ij}=
\begin{cases}
\widehat{H}_{ij},&i\neq j,\\
\dfrac{1}{2}+\widehat{H}_{ij},&i= j,
\end{cases}
\end{equation*}and
\begin{equation*}
\widehat{H}_{ij}=\int_{\partial\Omega_j}\dfrac{\partial w}{\partial n}\,\dd s_{j},\quad G_{ij}=\int_{\partial\Omega_j}w\,\dd s_{j},
\end{equation*}
we get 
\begin{equation*}
\sum_{j=1}^NH_{ij} u_j=\sum_{j=1}^N G_{ij}q_j.
\end{equation*}
Here we know some of $u_j$'s and some of $q_j$'s, i.e., we have exactly $N$
known quantities and $N$ unknowns. Then, by applying the boundary conditions to every midpoint of the elements, we have the equations in the matrix-vector form
\[HU=GQ.\]
Rearranging this system with known quantities on one side and the unknowns on the other, we get the linear system of equations $AX=B$ where $X$ is unknow values.
%%%%%%%%%%%%%%%%%%%%%
\begin{example}
Consider the Laplace equation in the rectangle plate with mixed boundary conditions as \cite{steady}
\begin{equation*}
 \begin{cases}
 \nabla^2 T=0,& 0<x<2 , 0<y<1,\\
 \dfrac{\partial u}{\partial n}(x, 0)=0, T(2,y)=0,\\
\dfrac{\partial u}{\partial n}(x, 1)=0,T(0,y)=1.
 \end{cases}
\end{equation*}
Using four constant element discretization with the collocation points at the center of the elements,  $H_{ij}$ and $G_{ij}$ for $i\neq j$ gives:
\begin{equation*}
\begin{split}
H_{ij} & =\int_{\partial\Omega_j}q^*\dd s_j =-\dfrac{1}{2\pi}\int _{X_j}^{X_{j+1}}-\dfrac{1}{ r^2}(r_xn_x+r_yn_y)\dd X\\
&=-\dfrac{l_e}{4\pi}\sum_{k=1}^4 \dfrac{w_k}{r_k^2}(r_{xk}n_{xk}+r_{yk}n_{yk}),\\
\end{split}
\end{equation*}
and
\begin{equation*}
\begin{split}
G_{ij}& =\int_{\partial\Omega_j}w\dd s_j  =\dfrac{1}{2\pi }\int_{X_j}^{X_{j+1}}\ln \left(\dfrac{1}{r}\right)\dd X\\
&=\dfrac{l_e}{4\pi}\sum_{k=1}^4 w_k\ln\left(\dfrac{1}{r_k}\right).
\end{split}
\end{equation*}
also for $i=j$, we have
$
H_{ij}=\dfrac{1}{2},
$
and
\begin{equation*}
G_{ij}=\dfrac{1}{\pi}\int_{0}^{\frac{l_e}{2}}\ln\left(\dfrac{1}{r}\right)\dd X=\dfrac{l_e}{2\pi}\left(\ln\left(\dfrac{2}{l_e}\right)+1\right),
\end{equation*}
where $l_e=X_{j+1}-X_j$ is the length of the element on which the integral is calculated, $w_k$ represent the weight functions from the Gauss integration scheme and $r_k$ is the distance from the $i^{th}$ source point to each integration points on the $j^{th}$ element. So
 \begin{small}
\begin{equation*}
H=
\begin{bmatrix}
0.5000& -0.1250& -0.2497 &-0.1250\\
-0.2113 &0.5000 &-0.2113& -0.0780\\
-0.2497& -0.1250 &0.5000 &-0.1250\\
-0.2113 &-0.0780& -0.2113& 0.5000
\end{bmatrix},
\end{equation*}
\begin{equation*}
G=
\begin{bmatrix}
0.3183& -0.0210& - 0.0421& -0.0210\\
-0.0176& 0.2695 &-0.0176 &- 0.1119\\
-0.0421& -0.0210 &0.3183 &-0.0210\\
-0.0176& -0.1119& -0.0176 &0.2695
\end{bmatrix}
\end{equation*}
\end{small}
By using from this equations unknow values in system $AX=B$ calculated.
\begin{equation*}
T_1=0.4993,\,T_3=0.4993,\, q_2=-0.7578,\, q_4=0.7575.
\end{equation*}
Consider for calculating temperature $T$ in interior points $$P_1=(0.75,0.5), P_2=(1,0.5), P_3=(1.5,0.5)$$   as
\begin{equation*}
T_{(i)}=\sum_{k=1}^NG_{ik}q_k -\sum_{k=1}^N H_{ik}T_k,
\end{equation*}
where
\begin{equation*}
G_{ik}=\int_{\Gamma _k}w_i\dd s_k, H_{ik}=\int_{\Gamma _k}\dfrac{\partial w_i}{\partial n}\dd s_k,
\end{equation*}
So
\begin{equation*}
T_{P_1}=0.5943, T_{P_2}=0.4869, T_{P_3}=0.3118,
\end{equation*}
and the exact values are
\[
T_{P_1}^e=0.625,T_{P_2}^e=0.5,T_{P_3}^e=0.25.
\]
Compared with the exact values , the solution shows quite a good approximation and stability for the internal points not only for the boundary points. Figure \ref{Fig1},\ref{Fig4}
\begin{figure}[h!]
\centering
\includegraphics[width=12cm]{laplac.jpg}
 \caption{The effect of increasing the number of elements on approximate solution Laplace equation}\label{Fig1}
\end{figure}
\begin{figure}[h!]
\centering
\includegraphics[width=12cm]{laplaczoom.eps}
\caption{The effect of increasing the number of elements on approximate solution Laplace equation (zoom)}\label{Fig4}
\end{figure}
\end{example}

%%%%%%%%%%%%%%%%
\section{Boundary element method for the helmholtz equation}
In this section we apply the boundary element method to the Helmholtz  equation \eqref{ma3} for $k\neq 0$. Consider weight function $w$ such that it satisfies
 \begin{equation}\label{helm}
  \nabla^2 w+k^2 w+\delta(X-X_0)=0,
 \end{equation}
 where $X_0=(\xi,\eta)$ is a singular point. The solutions  for the homogeneous case are the Hankel functions of the first and second kinds of order zero. So we write
\[w(r)=AH_0^{(2)}(kr),\]
where $r=\sqrt{(x-\xi)^2+(y-\eta)^2}$ and $H_0^{(2)}$  is the Hankel function of the second kind of the order zero and it  is representable in terms of the 0th-order Bessel functions of the first and second kinds $J_0$, $ Y_0$  as \[H_0^{(2)}(z)=J_0(z)-iY_0(z)=\bar{H}_0^{(1)}(z).\]
By using the small-argument approximation to the Hankel function as \cite{gibson}
\begin{equation*}
H_0^{(2)}(k r)\approx 1-i\dfrac{2}{\pi}\log(\dfrac{\gamma kr}{2}),\qquad\quad r\rightarrow 0,
\end{equation*}
where $\gamma=1.781$ and integrating \eqref{helm} over a very small circle of radius $\varepsilon$ centered at the origin, yielding
\begin{equation*}
A\int_\Omega [\nabla\cdot\nabla +k^2][1-i\dfrac{2}{\pi}\log(\dfrac{\gamma kr}{2})]\dd v =-1,
\end{equation*}
using the divergence theorem, the first term is converted to a line integral yielding
\begin{equation*}
-i\dfrac{2}{\pi}\int_0^{2\pi}\nabla\left(\log(\dfrac{\gamma k r}{2})\right)r\dd\phi=-4i,
\end{equation*}
the second term goes to zero, so
\[A=-\dfrac{i}{4},\]
and therefore we get the following fundamental solution for helmholtz equation as \cite{gibson}
\[w=\dfrac{1}{4i}H_0^{(2)}(kr),\]
Also we note that\cite{Regular}
\begin{equation*}
%\begin{split}
q^*=\dfrac{\partial w}{\partial n}=\dfrac{\partial w}{\partial r}\dfrac{\partial r}{\partial n}=-\dfrac{1}{4i}kH_1^{(2)}(kr)(\nabla r\cdot n)=-\dfrac{n_xr _x+n_y r _y}{4ir}kH_1^{(2)}(kr),
%\end{split}
\end{equation*}
where $r_x=x-\xi$ and $r_y=y-\eta$. To obtain the boundary integral equation, we integrate \eqref{ma3} over the domain $\Omega$ with a weighting function $ w$ to get,
\begin{equation*}
\int_\Omega (\nabla^2u +k^2u)w\dd v=0.
\end{equation*}
By using the second Green identity ,  gives
\begin{equation*}
\int_\Omega w\nabla^2 u\dd v=\int_\Omega u\nabla ^2w\dd v+\int_{\partial\Omega}(uq^*-wq)\dd s,
\end{equation*}
so
\begin{align}
\int_\Omega (\nabla^2 w+k^2w)u\dd v+\int_{\partial\Omega}(uq^*-wq)\dd s=0,
\end{align}
by using the properties of Dirac delta gives
\begin{equation*}
\int_\Omega (\nabla^2 w+k^2w)u\dd v =-u(X_0).
\end{equation*}
So the boundary integral equation for all the interior points is
\begin{equation*}
u(X_0)+\int_{\partial\Omega}u q^*\dd s=\int_{\partial\Omega}wq\dd s.
\end{equation*}
The boundary integral equation in general case is  \cite{wall}
\begin{equation*}
C(X_0)u(X_0)+\int_{\partial\Omega}u q^*\dd s=\int_{\partial\Omega}wq\dd s,
\end{equation*}
where $C(X_0)=\dfrac{\alpha}{2\pi}$ is a geometry-dependent parameter where $\alpha$  is internal angle. Assume that the  boundary is discretized and represented by $N$ linear constant elements as $\partial\Omega=\cup_{j=1}^N\partial\Omega_j$ and the source points situated on a smooth boundary ($\alpha =\pi$), the boundary integral equation   can be written as
\begin{equation*}
\dfrac{1}{2} u(X_0)+\int_{\cup_{j=1}^N\partial\Omega_j}u\dfrac{\partial w}{\partial n}\dd s=\int_{\cup_{j=1}^N\partial\Omega_j}w\dfrac{\partial u}{\partial n}\dd s.
\end{equation*}
Therefore
\begin{equation*}
\sum_{j=1}^NH_{ij} u_j=\sum_{j=1}^N G_{ij}q_j,
\end{equation*}
where
\begin{equation*}
H_{ij}=
\begin{cases}
\widehat{H}_{ij},&i\neq j,\\
\dfrac{1}{2}+\widehat{H}_{ij},&i= j,
\end{cases}
\end{equation*}
and
\begin{equation*}
\widehat{H}_{ij}=\int_{\partial\Omega_j}\dfrac{\partial w}{\partial n}\,\dd s_{j}, G_{ij}=\int_{\partial\Omega_j}w\,\dd s_{j}.
\end{equation*}
 By applying this boundary conditions to every midpoint of the elements, we have the equations in the matrix-vector form
$HU=GQ$ and rearranging this system, we get the linear system of equations $AX=B$ where $X$ is unknow values.
\begin{example}
Consider the Helmholtz equation \eqref{ma3} with $k=\dfrac{\sqrt{2}\pi}{4}$, the solution domain to be $0<x<1$, $0<y<1$, and the boundary conditions as \cite{Beginner}
\begin{equation*}
\begin{cases}
\nabla^2 u+k^2u=0, &0<x<1, 0<y<1,\\
u(0,y)=0,&0<y<1,\\
u(1,y)=\dfrac{\sqrt{2}}{2}\cos(\dfrac{\pi y}{4}),& 0<y<1,\\
\dfrac{\partial u}{\partial n}\bigg|_{y=0}=0,& 0<x<1,\\
\dfrac{\partial u}{\partial n}\bigg|_{y=1}=-\dfrac{\pi\sqrt{2}}{8}\sin(\dfrac{\pi x}{4}),& 0<x<1.
\end{cases}
\end{equation*}
The exact solution of this particular boundary value problem is
\begin{align*}
u(x,y)=\sin(\dfrac{\pi x}{4})\cos(\dfrac{\pi y}{4}).
\end{align*}
 Considering constant elements and applying the boundary element method for interior points, we get Figures \ref{Fig2}, \ref{Fig3}
 %\begin{table}
%\caption{  }
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline  \emph{Interior point}&\emph{Forty elements} &\emph{Eighty elements}&\emph{Exact solution}  \\
\hline $(0.1,0.1)$&$0.0879$&$0.0876$&$0.0782$ \\
\hline $(0.2,0.2)$&$0.1484$&$0.1490$ &$0.1545$    \\
\hline $(0.3,0.3)$& $0.2108$&$0.2121$&$0.2270$ \\
\hline $(0.4,0.4)$& $0.2733$& $0.2756$&$0.2939$  \\
\hline $(0.5,0.5)$& $0.3335$& $0.3375$&$0.3536$  \\
\hline $(0.6,0.6)$& $0.3880$& $0.3944$&$0.4045$ \\
\hline $(0.7,0.7)$& $0.4335$& $0.4418$&$0.4455$  \\
\hline $(0.8,0.8)$& $0.4678$& $0.4751$ &$0.4755$ \\
\hline $(0.9,0.9)$& $0.4900$& $0.4934$&$0.4938$  \\
\hline
\end{tabular}
\end{center}
%\end{table}
\begin{figure}[H]
\centering
\includegraphics[width=12cm]{jadidhelm1.jpg}
 \caption{The effect of increasing the number of elements on approximate solution helmholtz equation}\label{Fig2}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[width=9cm]{helmzooom.eps}
 \caption{The effect of increasing the number of elements on approximate solution helmholtz equation (zoom)}\label{Fig3}
\end{figure}
\end{example}

\bibliography{Ref1}
\bibliographystyle{amsplain}

 %\nocite{*}

%\bibliography{Ref2}
%\bibliographystyle{amsplain}
\end{document}