Asymptotic analysis of the Narrow Escape Problem in general shaped domain with several absorbing necks

When a Brownian particle is confined to a simply connected bounded domain $\Omega$ with a small absorbing window and the other part of the boundary is reflecting, it attempts to escape from the domain through this small absorbing window. The time that the Brownian particle takes from an initial position $x\in\Omega$ to escape via the absorbing window is called the mean first passage time(MFPT). The narrow escape problem(NEP) is to calculate the MFPT of the confined particle. Let $\Omega$ be a bounded simply connected domain in $\mathbb{R}^{2}$ with $C^{2}$-smooth boundary $\partial\Omega$. The boundary $\partial\Omega$ is composed of two disjoint parts, the absorbing part $\partial\Omega_{a}$ and the reflecting part $\partial\Omega_{r}$, such that $\partial\Omega=\partial\Omega_{a}\cup\partial\Omega_{r}$. We assume that $|\partial\Omega_{a}|=O(\epsilon)$ for $\epsilon\ll 1$, while $\partial\Omega$ is of order $1$. Suppose that a Brownian particle is confined at the position $x\in\Omega$. The MFPT $u(x)$, which depends on the starting position $x\in\Omega$, satisfies the following equation:

$$\begin{cases}\Delta u(x)=-\frac{1}{D}\quad&\mbox{in}~{}\Omega,\\[4.30554pt] \dfrac{\partial u}{\partial\nu}(x)=0&\mbox{on}~{}\partial\Omega_{r},\\[6.45831pt] u(x)=0&\mbox{on\ }\partial\Omega_{a},\end{cases}$$

(1.1)

where $\nu$ is the outer unit normal to $\partial\Omega$ and $D$ is the diffusion coefficient. In this paper, we consider the simplest form of pure diffusion with a unit diffusion coefficient $D=1$. The asymptotic analysis for NEP arises in deriving the asymptotic expansion of $u$ as $\epsilon\rightarrow 0$, from which one can estimate the MFPT of the confined particle. The above Dirichlet-Neumann type boundary value problem (1.1) can be derived by considering the probability density function of the Brownian particle at location $x$ at time $t$ through Fokker-Planck equation. More details on the derivation can be found in [2].

The NEP has attracted significant attention from the point view of mathematical and numerical concern due to its application in molecular biology and biophysics. There have been significant works in deriving the leading order term and higher order terms of the asymptotic expansion of MFPT in a regular bounded domain in two and three dimensions, see, for example, [2, 3, 7, 8, 10, 6, 16, 18]. Most existing results are focusing on the derivation of the MFPT in a domain with a single small absorbing gate. Whereas, a few works have been conducted for a domain with multiple small absorbing windows. The leading order term of the asymptotic formula of MFPT in a smooth bounded domain with multiple small absorbing windows has been studied in [1, 4, 5, 9, 11, 15].

In this paper, we rigorously derive the high leading order asymptotic expansion of the MFPT in a domain which is composed of a relatively big head and several narrow absorbing necks. In fact, the MFPT in such domain is closely related to the diffusion of particles in a cellular network since the cellular network are mostly composed of narrow necks connecting relatively larger head compartments. The shape and the distribution of the head compartments as well as the distribution of the tubules are involved in processes including active or diffusion transport of proteins, calcium signaling and so on [12]. How changes in shape of the compartment or the ratio of head to necks occur in response to specific cellular signals is important. However, how these structures regulate molecular trafficking and diffusion is unclear [17]. In [13], we considered a new model of Robin-Neumann type boundary value problem to compute the MFPT in a domain composed of a relatively big head and a single neck gate. Motivated by [13] and its possible implications on diffusion of particles in cellular networks, we consider the NEP problem in a domain which is composed of a relatively big head and several thin absorbing necks from mathematical point of view in this paper.

Let $\Omega_{h}$ be a simply connected bounded domain with a $C^{2}$ boundary. Let $\Omega_{n_{i}}$ be rectangular neck with length $L_{i}$ and width $2\epsilon_{i}$, where $\epsilon_{i}\ll 1$, for $i=1,\dots,N$. Each $\Omega_{n_{i}}$ is connected to $\Omega_{h}$ with connecting segment $\Gamma_{\epsilon_{i}}$. We assume that $\Gamma_{\epsilon_{i}}\cap\Gamma_{\epsilon_{j}}=\emptyset$, $i\neq j$, and moreover $dist(\Gamma_{\epsilon_{i}},\Gamma_{\epsilon_{j}})\geq c$, for some $c=O(1)$, in other words, they are well-separated. Let $\Omega=\Omega_{h}\cup\Omega_{n_{1}}\cup\dots\cup\Omega_{n_{N}}$. The boundary of $\Omega$ is divided into two parts, one is the reflecting boundary $\partial\Omega_{r}$, and the other absorbing boundary $\partial\Omega_{a}$, which consists of $N$ parts $\partial\Omega_{a_{i}}$, $i=1,,\dots,N$, where $\partial\Omega_{a_{i}}$ is the end of each neck. We assume that the size of each $\partial\Omega_{a_{i}}$ is much smaller than that of the whole boundary $\partial\Omega$. As an example, the geometric description for four necks is shown in Figure 1. The NEP is to calculate the MFPT $u$ which is the unique solution to the following Dirichlet-Neumann model:

$$\begin{cases}\Delta u(x)=-1\quad&\mbox{in}~{}\Omega,\\[4.30554pt] \dfrac{\partial u}{\partial\nu}(x)=0&\mbox{on}~{}\partial\Omega_{r},\\[6.45831pt] u(x)=0&\mbox{on\ }\partial\Omega_{a_{i}},\end{cases}$$

(1.2)

for all $i=1,\dots,N$, where $\nu$ is the outward unit normal to $\partial\Omega$. In this paper, we derive the high order asymptotic expansion of $u$ as $\max(\epsilon_{1},\dots,\epsilon_{N})\rightarrow 0$, from which one can estimate the escape time of the Brownian particle through the absorbing neck gate with high accuracy. In fact, three leading order term asymptotic solution to (1.2) is derived by means of solving an equivalent Neumann-Robin type boundary value problem, which is proposed in [13] to calculate the MFPT in a domain with a single absorbing neck window.

The main idea of the Neumann-Robin model is as follows: since the neck is relatively thin, it can be approximated as one dimensional, then the solution $u$ only varies along the neck direction. By dropping the neck and using the continuity of $u$ and its derivatives, it is technically equivalent as considering a Robin boundary condition. Following the same spirit, in this paper, we drop all the necks and instead consider Robin boundary condition at each connecting segment $\Gamma_{\epsilon_{i}}$, $i=1,\dots,N$. Hence, we can reformulate the original problem (1.2) as the following Neumann-Robin model:

$$\begin{cases}\Delta u(x)=-1\quad&\mbox{in}~{}\Omega_{h},\\[4.30554pt] \dfrac{\partial u}{\partial\nu}(x)=0&\mbox{on}~{}\partial\Omega_{r},\\[10.00002pt] \dfrac{\partial u}{\partial\nu}(x)+\alpha_{i}u(x)=\beta_{i}&\mbox{on}~{}\Gamma_{\epsilon_{i}}.\end{cases}$$

(1.3)

Here $\alpha_{i}$ and $\beta_{i}$ are constants which are determined in Section 3. We assume that each $\epsilon_{i}$ is sufficiently small such that $\alpha_{i}\epsilon_{i}\ll 1$.

By considering the Neumann-Robin model (1.3), we transform the problem from a singular domain $\Omega$ into a smooth domain $\Omega_{h}$ by dropping all the necks and assigning Robin condition on each connecting segment. In this way, we shall use layer potential techniques to derive the high order asymptotic solution $u$ to (1.3), and obtain the following main theorem of this paper.

Moreover, if $\Omega_{h}$ is a unit disk centered at $(0,0)$, the necks have the same length $L$ and the same width $2\epsilon$, then the MFPT $u$ has the following explicit formula up to order $O\left(\epsilon\ln^{2}\epsilon\right)$:

	$\displaystyle u\left(x\right)=$	$\displaystyle\frac{\left|\Omega_{h}\right|L}{2N}\frac{1}{\epsilon}-\frac{\left|\Omega_{h}\right|}{\pi N}\ln\epsilon-\frac{\left|\Omega_{h}\right|}{2\pi N}\left(2\ln 2-3\right)+\frac{L^{2}}{2}-\frac{2\left|\Omega_{h}\right|}{\pi N^{2}}\sum_{i=1}^{N-1}\sum_{j=i+1}^{N}\ln\left|s_{i}-s_{j}\right|$		(1.4)
		$\displaystyle+\frac{1}{4}\left(1-\left|x\right|^{2}\right)+\frac{\left|\Omega_{h}\right|}{\pi N}\sum_{i=1}^{N}{\ln\left|x-s_{i}\right|}+O\left(\epsilon\ln^{2}\epsilon\right),$

where $s_{i}$ is the center point on $\Gamma_{\epsilon_{i}}$. The first leading term of (1.4) is proportional to the length $L$ and of order $O(1/\epsilon)$, which is different from the well-known leading term $O(\ln\epsilon)$ for two dimensional NEP without necks. Due to the existence of narrow necks, it takes longer time for the particle to escape, which is natural and also interesting. Moreover, the three leading terms are explicitly given in formula (1.4). The formula (1.4) can approximate the MFPT with high accuracy, which can be seen from the numerical results in Section 5.

It is also worth mentioning that, let $N=1$, i.e., there is only a single neck gate, then the above formula becomes

	$\displaystyle u\left(x\right)$	$\displaystyle=\frac{\left|\Omega_{h}\right|L}{2}\frac{1}{\epsilon}-\frac{\left|\Omega_{h}\right|}{\pi}\ln\epsilon-\frac{\left|\Omega_{h}\right|}{2\pi}\left(2\ln 2-3\right)+\frac{L^{2}}{2}$		(1.5)
		$\displaystyle+\frac{1}{4}\left(1-\left|x\right|^{2}\right)+\frac{\left|\Omega_{h}\right|}{\pi}{\ln\left|x-s_{1}\right|}+O\left(\epsilon\right),$

which is exactly the same as the result in [13], where the MFPT was derived for a single neck exit.

By comparing (1.4) and (1.5), one can observe that the first and second leading order terms of (1.4), where there exist $N$ necks, is $1/N$ of that of (1.5), where there exists only a single neck, which is quite natural. But, the third leading order term $O(1)$ does not satisfy the $1/N$ relation between (1.4) and (1.5). The reason is that the $O(1)$ term depends not only on the location of the starting point, the length of the neck, but also on the location of the narrow necks as well as the nonlinear interaction between them. We will show how the $O(1)$ term contributes to the accuracy of the MFPT through numerical results in Section 5. As far as we know, this paper gives the first result on the three-term asymptotic expansion of the solution to NEP in a domain with several thin neck windows.

This paper is organized as follows. In section 2, we review the Neumann function for the Laplacian in $\mathbb{R}^{2}$. In section 3, we derive an equivalent Neumann-Robin model for the associated NEP in a domain with several thin necks. In section 4, we rigorously derive the high order asymptotic solution to the Neumann-Robin model by using layer potential techniques. Numerical experiments are provided in Section 5 to confirm the theoretical results. This paper ends with a short conclusion.

In this section, we review on the structure of Neumann function for a regular domain in $\mathbb{R}^{2}$ for further use, refer to [7].

Let $\Omega$ be a bounded domain in $\mathbb{R}^{2}$ with $C^{2}$ smooth boundary $\partial\Omega$, and let $G(x,z)$ be the Neumann function for $-\Delta$ in $\Omega$ with a given $z\in\Omega$. That is, $G(x,z)$ is the solution to the boundary value problem

$$\begin{cases}\Delta_{x}G(x,z)=-\delta_{z},&x\in\Omega,\\[4.30554pt] \displaystyle\frac{\partial G}{\partial\nu_{x}}=-\frac{1}{|\partial\Omega|},&x\in\partial\Omega,\\[10.00002pt] \displaystyle\int_{\partial\Omega}G(x,z)d\sigma(x)=0,&\end{cases}$$

where $\nu$ is the outer unit normal to the boundary $\partial\Omega$.

If $z\in\Omega$, then $G(x,z)$ can be written in the form

$\displaystyle G(x,z)=-\frac{1}{2\pi}\ln|x-z|+R_{\Omega}(x,z),\quad x\in\Omega,$

where $R_{\Omega}(x,z)$ is the regular part which belongs to $H^{3/2}(\Omega)$, the standard $L^{2}$ Sobolev space of order $3/2$, which solves the boundary value problem

$$\begin{cases}-\Delta_{x}R_{\partial\Omega}(x,z)=0,&x\in\Omega,\\[4.30554pt] \displaystyle\frac{\partial R_{\Omega}}{\partial\nu_{x}}\Big{|}_{x\in\partial\Omega}=-\frac{1}{|\partial\Omega|}+\frac{1}{2\pi}\frac{\langle x-z,\nu_{x}\rangle}{|x-z|^{2}},&x\in\partial\Omega,\end{cases}$$

where $\langle\cdot,\cdot\rangle$ denotes the inner product in $\mathbb{R}^{2}$.

If $z\in\partial\Omega$, then Neumann function on the boundary is denoted by $G_{\partial\Omega}$ and can be written as

$\displaystyle G_{\partial\Omega}(x,z)=-\frac{1}{\pi}\ln|x-z|+R_{\partial\Omega}(x,z),\quad x\in\Omega,~{}z\in\partial\Omega,$

(2.1)

where $R_{\partial\Omega}(x,z)$ has weaker singularity than $\ln|x-z|$ and solves the boundary value problem

$$\begin{cases}\Delta_{x}R_{\partial\Omega}(x,z)=0,&x\in\Omega,\\[4.30554pt] \displaystyle\frac{\partial R_{\partial\Omega}}{\partial\nu_{x}}\Big{|}_{x\in\partial\Omega}=-\frac{1}{|\partial\Omega|}+\frac{1}{\pi}\frac{\langle x-z,\nu_{x}\rangle}{|x-z|^{2}},&x\in\partial\Omega,~{}z\in\partial\Omega.\end{cases}$$

In particular, if $\Omega$ is the unit disk, then

$$\frac{\langle x-z,\nu_{x}\rangle}{|x-z|^{2}}=\frac{\langle x-z,x\rangle}{|x|^{2}+|z|^{2}-2x\cdot z}=\frac{1-x\cdot z}{2-2x\cdot z}=\frac{1}{2}$$

for $x\in\partial\Omega$. Therefore, $\frac{\partial R_{\partial\Omega}}{\partial\nu_{x}}\Big{|}_{x\in\partial\Omega}=0$, and hence $R_{\partial\Omega}(x,z)=$constant. Since $\int_{\partial\Omega}G(x,z)d\sigma(x)=0$, we have $R_{\partial\Omega}(x,z)=0$ for all $x\in\Omega$ and $z\in\partial\Omega$, and hence

$\displaystyle G_{\partial\Omega}(x,z)=-\frac{1}{\pi}\ln|x-z|,~{}x\in\Omega,~{}z\in\partial\Omega.$

We also have

$\displaystyle\int_{\Omega}G(x,z)dz=\frac{1}{4}(1-|x|^{2}).$

For later use, we introduce the integral operator $P:L^{2}[-1,1]\rightarrow L^{2}[-1,1]$, defined by

$$P[\phi](x)=\int_{-1}^{1}\ln|x-y|\phi(y)dy.$$

The operator $P$ is bounded (see Lemma 2.1 in [2]). The following quantity (obtained in [13]) is useful for later calculation

$$\int_{-1}^{1}{P\left[1\right]}\left(t\right)dt=4\ln 2-6.$$

(2.2)

In this section, we derive an equivalent Neumann-Robin type boundary value problem to (1.2) to solve the NEP in a composite domain which comprises a relatively large head and several narrow rectangular absorbing necks.

Let $\Omega_{h}$ be a simply connected bounded domain with a $C^{2}$ boundary. Let $\Omega_{n_{i}}$ be rectangular neck with length $L_{i}$ and width $2\epsilon_{i}$, where $\epsilon_{i}\ll 1$, for $i=1,\dots,N$. The necks are connected to the head domain $\Omega_{h}$ and they are well-separated. The connecting part between $\Omega_{h}$ and $\Omega_{n_{i}}$ is a small line segment which is denoted by $\Gamma_{\epsilon_{i}}$. Let $\partial\Omega_{a_{i}}$ be the exit, which is the end of the neck $\Omega_{n_{i}}$. Denote $\Omega:=\Omega_{h}\cup\Omega_{n_{1}}\cup\dots\cup\Omega_{n_{N}}$.

Consider the neck domain $\Omega_{n_{i}}$. Take the center point of $\Gamma_{\epsilon_{i}}$ as the point of origin $(0,0)$, then the neck direction is along the $x$- axis. Following (1.2), the MFPT $u$ in domain $\Omega_{n_{i}}$ should satisfy the boundary value problem:

$$\begin{cases}\Delta u(x)=-1\quad&\mbox{in}~{}\Omega_{n_{i}},\\[4.30554pt] \dfrac{\partial u}{\partial\nu}(x)=0&\mbox{on}~{}\partial\Omega_{r},\\[6.45831pt] u(x)=0&\mbox{on\ }\partial\Omega_{a_{i}},\\[2.15277pt] u(x)=u(x,y)&\mbox{on\ }\Gamma_{\epsilon_{i}},\end{cases}$$

where $x$ is the $x$- coordinate and $y$ is the $y$- coordinate in the two-dimensional domain.

As it is noted in paper [13], since the width of the neck is small, the neck domain $\Omega_{n_{i}}$ can be approximated as one-dimensional. The solution $u$ only varies along the $x$- coordinate, and thus solves the following ordinary differential equation:

		$\displaystyle\frac{d^{2}u(x,y)}{dx^{2}}=-1\quad~{}\mbox{for}~{}x\in[0,L_{i}],$
		$\displaystyle u(x,y)=0\quad\quad~{}\mbox{for}~{}x=0,$
		$\displaystyle u(x,y)=C_{i}\quad\quad~{}\mbox{for}~{}x=L_{i},$

where $C_{i}$ is constant. By separation of variables, we can solve the above ODE and obtain

$$u(x,y)=-\frac{1}{2}(L_{i}-x)^{2}+\left(\frac{C_{i}}{L_{i}}+\frac{L_{i}}{2}\right)(L_{i}-x),$$

(3.1)

where $x\in[0,L_{i}]$, $y\in[-\epsilon_{i},\epsilon_{i}]$.

Since $u$ and the derivative of $u$ are both continuous across the connecting segment $\Gamma_{\epsilon_{i}}$, they should satisfy a Robin type condition

$$\frac{\partial u}{\partial\nu}+\alpha_{i}u=\beta_{i}\quad\mbox{on}~{}\Gamma_{\epsilon_{i}},$$

(3.2)

where $\alpha_{i}$, $\beta_{i}$ are to be determined. In fact, substituting (3.1) into (3.2), one can see that

$$\alpha_{i}=\frac{1}{L_{i}},~{}\beta_{i}=\frac{L_{i}}{2}.$$

(3.3)

Therefore, by (3.2) and (3.3), the MFPT $u$ for a confined particle from the initial position $x$ in the head domain $\Omega_{h}$ to escape from the neck exits, is to solve the following equivalent Neumann-Robin boundary value problem:

$$\begin{cases}\Delta u(x)=-1\quad&\mbox{in}~{}\Omega_{h},\\[4.30554pt] \dfrac{\partial u}{\partial\nu}(x)=0&\mbox{on}~{}\partial\Omega_{r},\\[10.00002pt] \dfrac{\partial u}{\partial\nu}+\dfrac{1}{L_{i}}u=\dfrac{L_{i}}{2}&\mbox{on}~{}\Gamma_{\epsilon_{i}}.\end{cases}$$

(3.4)

Note that, by considering the Neumann-Robin model, we transform a singular domain $\Omega_{h}\cup\Omega_{n_{1}}\cup\dots\cup\Omega_{n_{N}}$ into a smooth domain $\Omega_{h}$ by dropping all the necks and assigning Robin condition on each connecting segment. In this way, we can derive high order asymptotic expansion of the MFPT by using layer potential techniques in the smooth domain.

In this section, we prove the main Theorem 1.1 in a smooth domain $\Omega_{h}$, where $\partial\Omega_{h}\in C^{2}(\mathbb{R}^{2})$. The high order asymptotic expansion of MFPT $u$ is rigorously derived by solving the equivalent Neumann-Robin model (3.4) using layer potential techniques.

Proof. Let

$$g(x)=\int_{\Omega_{h}}G(x,z)dz,~{}x\in\Omega_{h},$$

then $g$ solves the following boundary value problem

$$\begin{cases}\Delta g_{(}x)=-1\quad&\mbox{in}~{}\Omega_{h},\\[4.30554pt] \dfrac{\partial g}{\partial\nu}(x)=-\dfrac{\left|\Omega_{h}\right|}{\left|\partial\Omega_{h}\right|}&\mbox{on}~{}\partial\Omega_{h},\\[10.00002pt] \int_{\partial\Omega_{h}}{g}d\sigma=0.\end{cases}$$

(4.1)

Applying Green’s second formula, and using (3.4) and (4.1), we obtain

$$u\left(x\right)=g\left(x\right)+\sum_{i=1}^{N}{\int_{\Gamma_{\epsilon_{i}}}G_{\partial\Omega}\left(x,y\right){\frac{\partial u\left(y\right)}{\partial\nu}d\sigma\left(y\right)}}+C_{\epsilon},$$

(4.2)

where the constant $C_{\epsilon}$ is given by

$$C_{\epsilon}=\frac{1}{\left|\partial\Omega_{h}\right|}\int_{\partial\Omega_{h}}{u\left(y\right)d\sigma\left(y\right)}.$$

Here, $\frac{\partial u}{\partial\nu}$ on $\Gamma_{\epsilon_{i}}$, $i=1,2,\dots,N$, and $C_{\epsilon}$ are to be determined.

By (2.1), and the third and fourth Robin boundary conditions in (3.4), one can see that (4.2) can be written as

	$\displaystyle\frac{L_{j}^{2}}{2}-L_{j}\frac{\partial u\left(x\right)}{\partial\nu}=$	$\displaystyle g\left(x\right)-\frac{1}{\pi}\sum_{i=1}^{N}\int_{\Gamma_{\epsilon_{i}}}{\ln\left|x-y\right|\frac{\partial u\left(y\right)}{\partial\nu}d\sigma\left(y\right)}$		(4.3)
		$\displaystyle+\sum_{i=1}^{N}\int_{\Gamma_{\epsilon_{i}}}{R_{\partial\Omega_{h}}\left(x,y\right)\frac{\partial u\left(y\right)}{\partial\nu}d\sigma\left(y\right)}+C_{\epsilon},$

for $x\in\Gamma_{\epsilon_{j}}$, $j=1,2,\dots,N$, respectively.

Let $s_{i}$ be the center point on $\Gamma_{\epsilon_{i}}$, $i=1,2,\dots,N$, and let $x_{i}(t)$: $[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}]\rightarrow\mathbb{R}^{2}$ be the arc-length parametrization of $\Gamma_{\epsilon_{i}}$, i.e., $|x_{i}^{\prime}(t)|=1$ for all $t\in[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}]$. Then

$$\Gamma_{\epsilon_{i}}=\{x_{i}(t)|t\in[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}]\}.$$

For $i,j=1,2,\dots,N$, denote

		$\displaystyle f_{i}(t):=g(x_{i}(t)),~{}t\in[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}],$		(4.4)
		$\displaystyle r_{ij}(t,s):=R_{\partial\Omega_{h}}\left(x_{i}(t),x_{j}(s)\right),~{}t\in[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}],s\in[-\epsilon_{j}+s_{j},\epsilon_{j}+s_{j}],$
		$\displaystyle\phi_{i}(t):=\frac{\partial u(x_{i}(t))}{\partial\nu},~{}t\in[-\epsilon_{i}+s_{i},\epsilon_{i}+s_{i}].$

It then follows from (4.3) that

	$\displaystyle\frac{L_{j}^{2}}{2}-L_{j}\phi_{j}\left(t\right)$	$\displaystyle=f_{j}(t)-\frac{1}{\pi}\sum_{i=1}^{N}{\int_{s_{i}-\epsilon_{i}}^{s_{i}+\epsilon_{i}}{\ln\left|x_{j}\left(t\right)-x_{i}\left(s\right)\right|\phi_{i}\left(s\right)d\sigma\left(s\right)}}$		(4.5)
		$\displaystyle+\sum_{i=1}^{N}{\int_{s_{i}-\epsilon_{i}}^{s_{i}+\epsilon_{i}}r_{ji}(t,s)\phi_{i}\left(s\right)d\sigma\left(s\right)}+C_{\epsilon},$

where $t\in[-\epsilon_{j}+s_{j},\epsilon_{j}+s_{j}]$, $j=1,2,\dots,N$.

By changes of variables $t\rightarrow s_{i}+\epsilon_{i}t$ and $s\rightarrow s_{i}+\epsilon_{i}s$ for $t,s\in(s_{i}-\epsilon_{i},s_{i}+\epsilon_{i})$, the above integral equation (4.5) becomes

$$\left\{\begin{array}[]{l}\displaystyle\frac{L_{1}^{2}}{2}-\frac{L_{1}}{\epsilon_{1}}\phi_{1}\left(t\right)=f_{1}(s_{1}+\epsilon_{1}t)-\frac{1}{\pi}\sum\limits_{i=1}^{N}{\int_{-1}^{1}{\ln\left|x_{1}\left(s_{1}+\epsilon_{1}t\right)-x_{i}\left(s_{i}+\epsilon_{i}s\right)\right|\phi_{i}\left(s\right)ds}}\\ \displaystyle\hskip 113.81102pt+\sum\limits_{i=1}^{N}{\int_{-1}^{1}{r_{1i}(x_{1}\left(s_{1}+\epsilon_{1}t\right),x_{i}\left(s_{i}+\epsilon_{i}s\right))\phi_{i}\left(s\right)ds}}+C_{\epsilon},\\ \displaystyle\frac{L_{2}^{2}}{2}-\frac{L_{2}}{\epsilon_{2}}\phi_{2}\left(t\right)=f_{2}(s_{2}+\epsilon_{2}t)-\frac{1}{\pi}\sum\limits_{i=1}^{N}{\int_{-1}^{1}{\ln\left|x_{2}\left(s_{2}+\epsilon_{2}t\right)-x_{i}\left(s_{i}+\epsilon_{i}s\right)\right|\phi_{i}\left(s\right)ds}}\\ \displaystyle\hskip 113.81102pt+\sum\limits_{i=1}^{N}{\int_{-1}^{1}{r_{2i}(x_{2}\left(s_{2}+\epsilon_{2}t\right),x_{i}\left(s_{i}+\epsilon_{i}s\right))\phi_{i}\left(s\right)ds}}+C_{\epsilon},\\ \hskip 85.35826pt\vdots\\ \displaystyle\frac{L_{N}^{2}}{2}-\frac{L_{N}}{\epsilon_{N}}\phi_{N}\left(t\right)=f_{N}(s_{N}+\epsilon_{N}t)-\frac{1}{\pi}\sum\limits_{i=1}^{N}{\int_{-1}^{1}{\ln\left|x_{N}\left(s_{N}+\epsilon_{N}t\right)-x_{i}\left(s_{i}+\epsilon_{i}s\right)\right|\phi_{i}\left(s\right)ds}}\\ \displaystyle\hskip 113.81102pt+\sum\limits_{i=1}^{N}{\int_{-1}^{1}{r_{Ni}(x_{N}\left(s_{N}+\epsilon_{N}t\right),x_{i}\left(s_{i}+\epsilon_{i}s\right))\phi_{i}\left(s\right)ds}}+C_{\epsilon},\\ \end{array}\right.$$

(4.6)

where $\phi_{i}\left(t\right)=\epsilon_{i}\phi_{i}\left(s_{i}+\epsilon_{i}t\right)$ for $t\in[-1,1]$.

Define integral operators $P,P_{j}:L^{2}\left[-1,1\right]\rightarrow L^{2}\left[-1,1\right]$ as

	$\displaystyle P\left[\psi\right]$	$\displaystyle:=\int_{-1}^{1}{\ln\left|t-s\right|\psi\left(s\right)ds};$
	$\displaystyle P_{j}\left[\psi\right]$	$\displaystyle:=\frac{1}{\epsilon_{j}}\int_{-1}^{1}\big{\{}\ln\frac{\left|x_{j}\left(s_{j}+\epsilon_{j}t\right)-x_{j}\left(s_{j}+\epsilon_{j}s\right)\right|}{\epsilon_{j}\left|t-s\right|}+\pi r_{jj}\left(s_{j},s_{j}\right)$
		$\displaystyle\hskip 85.35826pt-\pi r_{jj}\left(s_{j}+\epsilon_{j}t,s_{j}+\epsilon_{j}s\right)\big{\}}\psi\left(s\right)ds,\quad j=1,\dots,N.$

One can easily check that

$$\left\{\begin{array}[]{l}\ln\left|x_{j}\left(s_{j}+\epsilon_{j}t\right)-x_{j}\left(s_{j}+\epsilon_{j}s\right)\right|=\ln\epsilon_{j}\left|t-s\right|+O\left(\epsilon_{j}\right),\\ \ln\left|x_{i}\left(s_{i}+\epsilon_{i}t\right)-x_{j}\left(s_{j}+\epsilon_{j}s\right)\right|=\ln\left|s_{i}-s_{j}\right|+O\left(\sqrt{\epsilon_{i}^{2}+\epsilon_{j}^{2}}\right),\\ \end{array}\right.$$

for $i,j=1,2,...,N;\;i\neq j$. Thus the operators $P$ and $P_{j}$ are bounded independently of $\epsilon_{j}$, $j=1,2,\dots,N$.

Integrating the first equation in (3.4) over $\Omega_{h}$ and using the divergence theorem, we obtain the compatibility condition

$$\sum_{i=1}^{N}\int_{\Gamma_{\epsilon_{i}}}{\partial_{\nu}ud\sigma}=-|\Omega_{h}|.$$

(4.7)

Let

$$C_{i}=\int_{-1}^{1}{\phi_{i}}\left(t\right)dt,\quad i=1,2,\dots,N.$$

(4.8)

By (4.7) and the third equation of (4.4), we have

$$\sum_{i=1}^{N}C_{i}=-|\Omega_{h}|.$$

(4.9)

Thus we can rewrite (4.6) as

$$\left\{\begin{array}[]{l}\left(I-\dfrac{\epsilon_{1}}{\pi L_{1}}\left(P+\epsilon_{1}P_{1}\right)\right)\phi_{1}\left(t\right)=\epsilon_{1}{C}_{\epsilon}^{1}+O\left(\epsilon_{1}\tilde{\epsilon}\right),\\[10.00002pt] \left(I-\dfrac{\epsilon_{2}}{\pi L_{2}}\left(P+\epsilon_{2}P_{2}\right)\right)\phi_{2}\left(t\right)=\epsilon_{2}{C}_{\epsilon}^{2}+O\left(\epsilon_{2}\tilde{\epsilon}\right),\\ \hskip 142.26378pt\vdots\\ \left(I-\dfrac{\epsilon_{N}}{\pi L_{N}}\left(P+\epsilon_{N}P_{N}\right)\right)\phi_{N}\left(t\right)=\epsilon_{N}{C}_{\epsilon}^{N}+O\left(\epsilon_{N}\tilde{\epsilon}\right),\end{array}\right.$$

(4.10)

where

$\displaystyle{C}_{\epsilon}^{i}=\frac{L_{i}}{2}-\frac{1}{L_{i}}\left(f_{i}\left(s_{i}\right)+\sum_{j=1}^{N}C_{j}r_{ij}+C_{\epsilon}\right)+\frac{1}{\pi L_{i}}\sum_{j\neq i}C_{j}\ln\left|s_{i}-s_{j}\right|+\frac{C_{i}}{\pi L_{i}}\ln\epsilon_{i},$

(4.11)

and $\tilde{\epsilon}=\max\{\epsilon_{1},\epsilon_{2},\dots,\epsilon_{N}\}$, $r_{ij}=r_{ij}(s_{i},s_{j})$ for $i,j=1,2,\dots,N$. Note that $r_{ij}=r_{ji}$, when $i\neq j$.

Since we assume that $\frac{\epsilon_{i}}{\pi L_{i}}\ll 1$ for $i=1,\dots,N$, one can easily see that

$$\left(I-\frac{\epsilon_{i}}{\pi L_{i}}\left(P+\epsilon_{i}P_{i}\right)\right)^{-1}=I+\frac{\epsilon_{i}}{\pi L_{i}}\left(P+\epsilon_{i}P_{i}\right)+O\left(\epsilon_{i}^{2}\right).$$

(4.12)

Thus, by (4.10) and (4.12), we have

$$\phi_{i}\left(t\right)=\epsilon_{i}{C}_{\epsilon}^{i}\left(1+\dfrac{\epsilon_{i}}{\pi L_{i}}P\left[1\right]\left(t\right)\right)+O\left(\epsilon_{i}\tilde{\epsilon}\right).$$

(4.13)

Integrating (4.13) over $(-1,1)$ and by (4.8), we have

$$\dfrac{C_{i}}{2\epsilon_{i}}={C}_{\epsilon}^{i}\left(1+\dfrac{\epsilon_{i}}{2\pi L_{i}}\int_{-1}^{1}{P\left[1\right]}\left(t\right)dt\right)+O\left(\tilde{\epsilon}\right).$$

(4.14)

It is easy to see that

$$\left(1+\frac{\epsilon_{i}}{2\pi L_{i}}\int_{-1}^{1}{P\left[1\right]}\left(t\right)dt\right)^{-1}=1-\frac{\epsilon_{i}}{2\pi L_{i}}\int_{-1}^{1}{P\left[1\right]}\left(t\right)dt+O\left(\epsilon^{2}_{i}\right).$$

Hence, by (2.2) and (4.14) we have

$${C}_{\epsilon}^{i}=\frac{C_{i}}{2\epsilon_{i}}-\frac{C_{i}}{2\pi L_{i}}\left(2\ln 2-3\right)+O\left(\tilde{\epsilon}\right),\quad i=1,2,\dots,N.$$

(4.15)

Comparing (4.15) with (4.11), together with (4.9), we obtain the system of $C_{1},\dots,C_{N}$ and $C_{\epsilon}$:

$$\mathcal{K}\left[\begin{array}[]{c}C_{1}\\ \vdots\\ C_{N}\\ C_{\epsilon}\end{array}\right]=\left[\begin{array}[]{c}\frac{L_{1}^{2}}{2}+f_{1}\left(s_{1}\right)+O\left(\tilde{\epsilon}\right)\\ \vdots\\ \frac{L_{N}^{2}}{2}+f_{N}\left(s_{N}\right)+O\left(\tilde{\epsilon}\right)\\ -\left|\Omega_{h}\right|\\ \end{array}\right],$$

(4.16)

where

$$\mathcal{K}=\left[\begin{matrix}A_{\epsilon}^{1}&\frac{-\ln\left|s_{1}-s_{2}\right|}{\pi}+r_{12}&\frac{-\ln\left|s_{1}-s_{3}\right|}{\pi}+r_{13}&\cdots&\frac{-\ln\left|s_{1}-s_{N}\right|}{\pi}+r_{1N}&1\\ \frac{-\ln\left|s_{2}-s_{1}\right|}{\pi}+r_{12}&A_{\epsilon}^{2}&\frac{-\ln\left|s_{2}-s_{3}\right|}{\pi}+r_{23}&...&\frac{-\ln\left|s_{2}-s_{N}\right|}{\pi}+r_{2N}&1\\ \frac{-\ln\left|s_{3}-s_{1}\right|}{\pi}+r_{13}&\frac{-\ln\left|s_{3}-s_{2}\right|}{\pi}+r_{23}&A_{\epsilon}^{3}&\cdots&\frac{-\ln\left|s_{3}-s_{N}\right|}{\pi}+r_{3N}&1\\ \vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ \frac{-\ln\left|s_{N}-s_{1}\right|}{\pi}+r_{1N}&\frac{-\ln\left|s_{N}-s_{2}\right|}{\pi}+r_{2N}&\frac{-\ln\left|s_{N}-s_{3}\right|}{\pi}+r_{3N}&\cdots&A_{\epsilon}^{N}&1\\ 1&1&1&\cdots&1&0\\ \end{matrix}\right]$$

and $A_{\epsilon}^{i}=\frac{L_{i}}{2\epsilon_{i}}-\frac{\ln\epsilon_{i}}{\pi}-\frac{2\ln 2-3}{2\pi}+r_{ii}$. Solving (4.16), we obtain

$$C_{i}=\frac{-\left|\Omega_{h}\right|\epsilon_{i}/L_{i}}{\sum\limits_{k=1}^{N}{\epsilon_{k}/L_{k}}}+O(\tilde{\epsilon}\ln\epsilon_{i}),$$

and

$$C_{\epsilon}=\frac{\left|\Omega_{h}\right|}{2\sum_{i=1}^{N}{\epsilon_{i}/L_{i}}}+\frac{\left|\Omega_{h}\right|}{\pi}\left(\sum_{i=1}^{N-1}{\sum_{j=i+1}^{N}{T_{ij}\ln\epsilon_{i}\epsilon_{j}}}-\sum_{i=1}^{N}{F_{i}\ln\epsilon_{i}}\right)+O\left(1\right),$$

(4.17)

where the constants $T_{ij}$, $F_{i}$ are bounded independently of $\epsilon_{i}$ that are given by

$$T_{ij}=\frac{\epsilon_{i}\epsilon_{j}}{L_{i}L_{j}}\Big{/}\left(\sum_{k=1}^{N}{\frac{\epsilon_{k}}{L_{k}}}\right)^{2},\quad F_{i}=\frac{\epsilon_{i}}{L_{i}}\Big{/}\sum_{k=1}^{N}{\frac{\epsilon_{k}}{L_{k}}}.$$

(4.18)

Thus, by (4.13), (4.15) and the third equation of (4.4), we have

$$\frac{\partial u\left(x_{i}(t)\right)}{\partial\nu}=\frac{C_{i}}{2\epsilon_{i}}-\frac{C_{i}}{2\pi L_{i}}\left(2\ln 2-3\right)+\frac{C_{i}}{2\pi L_{i}}P\left[1\right]\left(\frac{t-s_{i}}{\epsilon_{i}}\right)+O\left(\tilde{\epsilon}\right),$$

(4.19)

where $t\in\left(s_{i}-\epsilon_{i},s_{i}+\epsilon_{i}\right),\;i=1,\dots,N$ .