Numerical computation of the capacity of generalized condensers

TL;DR

This paper introduces a boundary integral method for numerically calculating the capacity of generalized condensers with complex geometries, also enabling harmonic measure computation in multiply connected domains.

Contribution

The paper presents a versatile boundary integral approach applicable to various complex geometries of generalized condensers, including piecewise smooth and rectilinear configurations.

Findings

Effective numerical computation of condenser capacity for complex geometries

Method applicable to multiply connected domains for harmonic measure calculation

Broad applicability to different boundary shapes

Abstract

We present a boundary integral method for numerical computation of the capacity of generalized condensers. The presented method applies to a wide variety of generalized condenser geometry including the cases when the plates of the generalized condenser are bordered by piecewise smooth Jordan curves or are rectilinear slits. The presented method is used also to compute the harmonic measure in multiply connected domains.

Equations

$$\mathrm{cap}(C)=\underset{G}{\iint }|\nabla u{|}^{2}dxdy$$

$$\Gamma =\partial G=\bigcup _{k=1}^{m+\ell }{\Gamma }_k,$$

$$m^{\prime}=\begin{cases}m-1,&\mbox{if $G$ is bounded and $B$ is unbounded},\\ m,&\mbox{otherwise},\end{cases}$$

$$\ell^{\prime}=\begin{cases}\ell-1,&\mbox{if $G$ is bounded and $B$ is bounded},\\ \ell,&\mbox{otherwise}.\end{cases}$$

$$u(\zeta )$$

$$\frac{\partial u}{\partial \mathbf{n}}(\zeta )$$

$$F(z)=g(z)-\sum _{k=1}^{m^{\prime }}{a}_k\log (z-{\alpha }_k)-\sum _{k=1}^{{\ell }^{\prime }}{b}_k\log (z-{\beta }_k)$$

$$a_k=\frac{1}{2\pi }{\int }_{{\Gamma }_k}\frac{\partial u}{\partial \mathbf{n}}ds,k=1,2,\dots ,m^{\prime },$$

$$b_k=\frac{1}{2\pi }{\int }_{{\Gamma }_{m+k}}\frac{\partial u}{\partial \mathbf{n}}ds,k=1,2,\dots ,{\ell }^{\prime }.$$

$$F(z)=g(z)-\sum _{k=1}^{m^{\prime }}{a}_k\log (z-{\alpha }_k).$$

$${\int }_{\Gamma }\frac{\partial u}{\partial \mathbf{n}}ds=0,$$

$$\sum _{k=1}^m{\int }_{{\Gamma }_k}\frac{\partial u}{\partial \mathbf{n}}ds=0.$$

$$a_m=\frac{1}{2\pi }{\int }_{{\Gamma }_m}\frac{\partial u}{\partial \mathbf{n}}ds.$$

$$\sum _{k=1}^m{a}_k=\sum _{k=1}^m\frac{1}{2\pi }{\int }_{{\Gamma }_k}\frac{\partial u}{\partial \mathbf{n}}ds=0,$$

$$a_m=-\sum _{k=1}^{m-1}{a}_k.$$

$$\mathrm{cap}(C)={\int }_{\partial G}u\frac{\partial u}{\partial \mathbf{n}}ds.$$

$$\mathrm{cap}(C)=\sum _{k=1}^m{\delta }_k{\int }_{{\Gamma }_k}\frac{\partial u}{\partial \mathbf{n}}ds=2\pi \sum _{k=1}^m{\delta }_k{a}_k.$$

$$J=\underset{k=1}{\overset{m+\ell }{⨆}}{J}_k=\bigcup _{k=1}^{m+\ell }\{(t,k):t\in J_k\}.$$

$$\eta (t,k)={\eta }_k(t),t\in J_k,k=1,2,\dots ,m+\ell .$$

$$\eta(t)=\left\{\begin{array}[]{l@{\hspace{0.5cm}}l}\eta_{1}(t),\hfil\hskip 14.22636pt&t\in J_{1},\\ \eta_{2}(t),\hfil\hskip 14.22636pt&t\in J_{2},\\ \hskip 8.5359pt\vdots\hfil\hskip 14.22636pt\\ \eta_{m+\ell}(t),\hfil\hskip 14.22636pt&t\in J_{m+\ell}.\end{array}\right.$$

$$\mathrm{Re}[F(\eta (t))]={\delta }_k,\eta (t)\in {\Gamma }_k,k=1,2,\dots ,m.$$

$$\mathbf{n}(\eta (t))=-\mathrm{i}\mathbf{T}(\eta (t))=-\mathrm{i}\frac{{\eta }^{\prime }(t)}{|{\eta }^{\prime }(t)|}=e^{\mathrm{i}\nu (\eta (t))}$$

$$\frac{\partial u}{\partial \mathbf{n}}=\nabla u\cdot \mathbf{n}=\cos (\nu )\frac{\partial u}{\partial x}+\sin (\nu )\frac{\partial u}{\partial y}=\mathrm{Re}\left[e^{\mathrm{i}\nu }\left(\frac{\partial u}{\partial x}-\mathrm{i}\frac{\partial u}{\partial y}\right)\right]=\mathrm{Re}\left[\frac{-\mathrm{i}{\eta }^{\prime }(t)}{|{\eta }^{\prime }(t)|}{F}^{\prime }(\eta (t))\right]$$

$$\mathrm{Re}\left[-\mathrm{i}{\eta }^{\prime }(t)F^{\prime }(\eta (t))\right]=0,\eta (t)\in {\Gamma }_{m+k},k=1,2,\dots ,\ell .$$

$$\mathrm{Re}[-\mathrm{i}F(\eta (t))]={\nu }_k,\eta (t)\in {\Gamma }_{m+k},k=1,2,\dots ,\ell ,$$

$$\mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}F(\eta (t))\right]=\delta (t)+\nu (t)$$

$$\theta (t)=\{\begin{array}{ll}0,& t\in J_1,\\ ⋮& \\ 0,& t\in J_m,\\ \pi /2,& t\in J_{m+1},\\ ⋮& \\ \pi /2,& t\in J_{m+\ell },\end{array},\delta (t)=\{\begin{array}{ll}{\delta }_1,& t\in J_1,\\ ⋮& \\ {\delta }_m,& t\in J_m,\\ 0,& t\in J_{m+1},\\ ⋮& \\ 0,& t\in J_{m+\ell },\end{array},\nu (t)=\{\begin{array}{ll}0,& t\in J_1,\\ ⋮& \\ 0,& t\in J_m,\\ {\nu }_1,& t\in J_{m+1},\\ ⋮& \\ {\nu }_{\ell },& t\in J_{m+\ell },\end{array}$$

$$\mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}g(\eta (t))\right]=\delta (t)+\nu (t)+\sum _{k=1}^{m^{\prime }}{a}_k\mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}\log (\eta (t)-{\alpha }_k)\right].$$

$${\gamma }_k(t)=\{\begin{array}{ll}\mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}\log (\eta (t)-{\alpha }_k)\right],& \mathrm{if}{\ell }^{\prime }=\ell ,\\ \mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}\log \frac{\eta (t)-{\alpha }_k}{\eta (t)-\alpha }\right],& \mathrm{if}{\ell }^{\prime }=\ell -1,\end{array}$$

$$\sum _{k=1}^{m^{\prime }}{a}_k{\gamma }_k(t)=\sum _{k=1}^{m^{\prime }}{a}_k\mathrm{Re}\left[e^{-\mathrm{i}\theta (t)}\log (\eta (t)-{\alpha }_k)\right].$$

$$\sum _{k=1}^{m^{\prime }}{a}_k{\gamma }_k(t)$$

Full text

††footnotetext: File: gc-2019-8-6.tex, printed: 2024-3-2, 15.29

1. Introduction

The conformal capacity of condensers is an important notion in geometric function theory [ah, avv, dek, du, kuh, ps, Vas02] and in various applications of electronics. However, the analytic forms of the capacity are known only for special types of condensers. So, the use of numerical methods for computing such capacity is unavoidable. Indeed, numerical computing of capacity of condensers have been intensively studied in the literature, see e.g., [bbgghv, bsv, hrv1, hrv2, dek]. The capacity of condensers is one of the several “conformal invariants” which are powerful tools in complex analysis. Some of the other important examples of conformal invariants are the harmonic measure, the logarithmic capacity, the extremal length, the reduced extremal length, and the hyperbolic distance [ah, du, du2, gm, Vas02, vuo]. Numerical computing of such invariants has been studied also in the literature, see e.g., [bbgghv, LSN17, ps, Ran2, RR].

Capacity of generalized condensers is another important example of conformal invariants [du, de, dk1, dk2, Vas02]. In this paper, we present a numerical method for computing the capacity of generalized condensers. We consider the case in which the plates of the generalized condensers are bordered by piecewise smooth Jordan curves or are rectilinear slits. As far as we know, the proposed method is the first numerical method for computing the numerical values of the capacity of the generalized condensers. The boundary integral equation with the generalized Neumann kernel [Nas-ETNA, Weg-Nas] plays a key role in developing our method. The presented method can be used also to compute the harmonic measure in multiply connected domains.

Let $B$ be an open subset of $\overline{{\mathbb{C}}}={\mathbb{C}}\cup\{\infty\}$. We consider generalized condensers of the form $C=(B,E,\delta)$ where $E=\{E_{k}\}_{k=1}^{m}$, $m\geq 2$, is a collection of nonempty closed pairwise disjoint sets in $E_{k}\subset B$ and $\delta=\{\delta_{k}\}_{k=1}^{m}$ is a collection of real numbers containing at least two different numbers. The set $G=B\backslash\cup_{k=1}^{m}E_{k}$ is called the field of the condenser $C$, the sets $E_{k}$ are the plates of the condenser, and the $\delta_{k}$ are the levels of the potential of the plates $E_{k}$, $k=1,2,\ldots,m$ [du, p. 12]. We assume that $G$ is a finitely connected domain without isolated boundary points and that $\partial G\cap\left(\cup_{k=1}^{m}E_{k}\right)$ consists of $m$ piecewise smooth Jordan curves, then the conformal capacity of $C$, $\mathrm{cap}(C)$, is given by the Dirichlet integral [du, p. 13, p. 305]

[TABLE]

where $u$ is the potential function of the condenser $C$, i.e., $u$ is continuous in $\overline{G}$, harmonic in $G$, and equal to $\delta_{k}$ on $\partial E_{k}$ for $k=1,2,\ldots,m$ and satisfies $\partial u/\partial{\bf n}=0$ on $\partial B\backslash\cup_{k=1}^{m}E_{k}$ where $\partial u/\partial{\bf n}$ denotes the directional derivative of $u$ along the outward normal.

The analytical description of the problem is given in Section 2 and it is based on the classical theory of integral equations [Mik64] and on the definition of the generalized capacity due to Dubinin [du]. In Section 3 we formulate the computational problem as a Riemann-Hilbert problem and prove a preliminary analytical result. The main theoretical results are presented in Section 4 and they deal with unique solvability of algebraic linear systems related to the Riemann-Hilbert problem. Also an outline of an algorithm for the numerical solution of the integral equation is given. In Section 5 we give a MATLAB implementation of the algorithm. This algorithm is tested in Section 6 in the case of capacity computation of condensers with piecewise smooth boundary curves and results are compared, with good agreement of results, to earlier numerical results from [bsv]. In Section 7 we apply the algorithm for the computation of the capacity of generalized condensers. In Section 8, we use the presented algorithm with the help of conformal mappings to compute the capacity of rectilinear slit condensers. In the final Section 9 we show that the same method also works for the computation of the harmonic measure.

2. The potential function

In this paper, for $k=1,2,\ldots,m$, we assume that $E_{k}=\overline{G_{k}}$ where $G_{k}$ is a simply connected domain bordered by a piecewise smooth Jordan curve $\Gamma_{k}$. We assume also that either $B={\mathbb{C}}$ or $B\subsetneq{\mathbb{C}}$ is a multiply connected domain of connectivity $\ell\geq 1$ bordered by $\ell$ piecewise smooth Jordan curves $\Gamma_{k}$ for $k=m+1,m+2,\ldots,m+\ell$. We assume $\ell=0$ when $B={\mathbb{C}}$ and $\infty\in B$ when $B$ is unbounded. Then, the field of the condenser is the multiply connected domain $G$ of connectivity $m+\ell$ bordered by

[TABLE]

where the orientation of the curves $\Gamma_{k}$ is such that $G$ is always on the left of $\Gamma_{k}$ for $k=1,2,\ldots,m+\ell$. For each $k=m+1,m+2,\ldots,m+\ell$, the simply connected domain on the right of $\Gamma_{k}$ will be denoted by $G_{k}$.

The domain $G$ is either bounded or unbounded. If $G$ is unbounded, we assume $\infty\in G$. If $G$ is bounded, then one of the simply connected domains $G_{1},\ldots,G_{m}$ or $G_{m+1},\ldots,G_{m+\ell}$ is unbounded and contains $\infty$. If the unbounded domain is one of the domains $G_{1},\ldots,G_{m}$, then we assume it is the domain $G_{m}$ and the curve $\Gamma_{m}$ enclose all the other curves $\Gamma_{1},\ldots,\Gamma_{m-1},\Gamma_{m+1},\ldots,\Gamma_{m+\ell}$. Similarly, if the unbounded domain is one of the domains $G_{m+1},\ldots,G_{m+\ell}$, then we assume it is the domain $G_{m+\ell}$ and the curve $\Gamma_{m+\ell}$ enclose all the other curves $\Gamma_{1},\ldots,\Gamma_{m+\ell-1}$. Based on the boundedness of the domains $B$ and $G$, we define the integers $m^{\prime}$ and $l^{\prime}$ by

[TABLE]

and

[TABLE]

In particular, if $G$ is unbounded, then $B$ is unbounded (since $G\subseteq B$), $m^{\prime}=m$, $\ell^{\prime}=\ell$, and hence $m^{\prime}+\ell^{\prime}=m+\ell$. If $G$ is bounded, then either $m^{\prime}=m-1$ or $\ell^{\prime}=\ell-1$ and hence $m^{\prime}+\ell^{\prime}=m+\ell-1$. Further, $m^{\prime}=m-1$ means that $\Gamma_{m}$ is the external boundary component of $G$. Similarly, $\ell^{\prime}=\ell-1$ means that the external boundary component of $G$ is $\Gamma_{m+\ell}$. With these definitions of $m^{\prime}$ and $\ell^{\prime}$, the domains $G_{1},\ldots,G_{m^{\prime}}$ and $G_{m+1},\ldots,G_{m+\ell^{\prime}}$ are bounded simply connected domains. For each of these bounded domains, we assume that $\alpha_{k}$ is an auxiliary point in $G_{k}$ for each $k=1,2,\ldots,m^{\prime}$ and $\beta_{k}$ is an auxiliary point in $G_{m+k}$ for each $k=1,2,\ldots,\ell^{\prime}$.

The potential function $u$ is then a solution of the Laplace equation $\Delta u=0$ with the mixed Dirichlet-Neumann boundary condition

[TABLE]

Note that the boundary value problem (2.3) reduces to a Dirichlet problem for $\ell=0$. Note also that the problem (2.3) does not reduce to a Neumann problem since $m\geq 2$. The problem (2.3) has a unique solution $u$ [IS19].

A more general form of such mixed boundary value problem has been considered in [IS19] using a Cauchy integral method and in [Nas-bvp, Nas-jam] using the boundary integral equation with the generalized Neumann kernel. Due to the simple forms of the boundary conditions in (2.3), the method presented in [Nas-bvp, Nas-jam] will be further simplified in this paper to obtain a simple, fast, and accurate method for computing the potential function $u$ and the capacity $\mathrm{cap}(C)$ of the generalized condenser $C$.

The harmonic function $u$ is the real part of an analytic function $F$ in $G$. The function $F$ is not necessarily single-valued, but it can be written as [Gak, gm, Mik64, Mus]

[TABLE]

where $g$ is a single-valued analytic function in $G$ and $a_{1},\ldots,a_{m^{\prime}},b_{1},\ldots,b_{\ell^{\prime}}$ are undetermined real constants such that [Mik64, §31]

[TABLE]

and

[TABLE]

Hence, using (2.3b), we have $b_{k}=0$ for all $k=1,2,\ldots,\ell^{\prime}$. Thus, the function $F$ has the representation

[TABLE]

Since $u$ is harmonic in the domain $G$, then [Mik64]

[TABLE]

which in view of (2.3b) implies that

[TABLE]

Recall that $a_{1},\ldots,a_{m^{\prime}}$ are given in (2.5). So, if $m^{\prime}=m-1$, we define

[TABLE]

Hence, it follows from (2.5), (2.7), and (2.8) that

[TABLE]

which implies, in the case $m^{\prime}=m-1$, that

[TABLE]

Using Green’s formula [du, p. 4], Equation (1.1) can be written as

[TABLE]

Since $\partial u/\partial{\bf n}=0$ on $\partial B=\cup_{k=1}^{\ell}\Gamma_{m+k}$ and $u=\delta_{k}$ on $\Gamma_{k}$ for $k=1,2,\ldots,m$, then in view of (2.5) and (2.8), we have

[TABLE]

Equation (2.12) gives us a simple formula for computing the capacity of the generalized condenser $C$ in terms of the levels $\delta_{k}$ of the potential of the plates and the values of the constants $a_{k}$ for $k=1,2,\ldots,m$.

In this paper, the boundary integral equation with the generalized Neumann kernel will be used to compute the constants $a_{k}$ as well as the values of the function $u(z)$ for $z\in G$. However, to use the integral equation, we will first reformulate the above mixed boundary value problem as a Riemann-Hilbert problem as it will be described in the next section. Solving the mixed boundary value problem by reducing it to a Riemann-Hilbert problem is a well known approach and has been used by many researchers in the literature (see e.g., [Nas-bvp, Gak, Haas, Mus, Nas-jam]).

3. The Riemann-Hilbert problem

For each $k=1,2,\ldots,m+\ell$, the boundary component $\Gamma_{k}$ is parametrized by a $2\pi$-periodic complex function $\eta_{k}(t)$, $t\in J_{k}:=[0,2\pi]$. The total parameter domain $J$ is the disjoint union of the $m+\ell$ intervals $J_{1},\ldots,J_{m+\ell}$,

[TABLE]

The elements of $J$ are ordered pairs $(t,k)$ where $k$ is an auxiliary index indicating which of the intervals contains the point $t$ [Nas-ETNA]. A parametrization of the whole boundary $\Gamma$ is then defined by

[TABLE]

For a given $t$, the value of auxiliary index $k$ such that $t\in J_{k}$ will be always clear from the context. So we replace the pair $(t,k)$ in the left-hand side of (3.1) by $t$ in the same way as in [Nas-ETNA]. Thus, the function $\eta$ in (3.1) is written as

[TABLE]

Since $u=\delta_{k}$ is known on the boundary components $\Gamma_{k}$ for $k=1,2,\ldots,m$ and since $u={\rm Re}\,F$, then the boundary values of the function $F$ satisfy

[TABLE]

On the boundaries $\partial B=\cup_{k=1}^{\ell}\Gamma_{m+k}$, the potential function $u$ satisfies the boundary condition $\partial u/\partial{\bf n}=0$ where ${\bf n}$ is the outward normal vector on $\partial B$. Let ${\bf T}$ be the unit tangent vector on $\partial B$. Then, for $\eta(t)\in\partial B$,

[TABLE]

where $\nu(\eta(t))$ is the angle between the positive real axis and the normal vector ${\bf n}(\eta(t))$. Using the Cauchy-Riemann equations, the derivative of the analytic function $F$ is then $F^{\prime}(z)=\frac{\partial u(z)}{\partial x}-\mathrm{i}\frac{\partial u(z)}{\partial y}$. Thus,

[TABLE]

which, in view of (2.3b), implies that

[TABLE]

Integrating with respect to the parameter $t$ yields

[TABLE]

where $\nu_{1},\nu_{2},\ldots,\nu_{\ell}$ are real constants of integration. Thus, by (3.3) and (3.6), the boundary values of the function $F$ satisfy the boundary condition

[TABLE]

where

[TABLE]

i.e., $\theta(t)=0$ and $\nu(t)=0$ for $\ell=0$. Then, it follows from (2.6) that the single-valued analytic function $g$ satisfies the boundary condition

[TABLE]

Lemma 3.9**.**

The functions $\gamma_{k}$, for $k=1,\ldots,m^{\prime}$, defined on $J$ by

[TABLE]

are periodic for $t\in J_{j}$, $j=1,2,\ldots,m+\ell$. For both cases, we have

[TABLE]

Proof.

Since $\theta(t)=0$ when $t\in J_{j}$ for each $j=1,2,\ldots,m$, then the functions $\gamma_{k}(t)$ in (3.10) are periodic for $t\in J_{j}$ for each $j=1,2,\ldots,m$.

When $t\in J_{j}$ for each $j=m+1,m+2,\ldots,m+\ell$, we have the following two cases:

a) $\ell^{\prime}=\ell$. For this case, $\Gamma_{m+\ell}$ is not the external boundary component of $G$. Recall that, for each $k=1,2,\ldots,m^{\prime}$, $\alpha_{k}$ is in the interior of the curve $\Gamma_{k}$. Thus, none of the auxiliary points $\alpha_{1},\ldots,\alpha_{m^{\prime}}$ is interior to any of the curves $\Gamma_{m+1},\ldots,\Gamma_{m+\ell}$. Hence, the winding number of the function $z-\alpha_{k}$ is always zero along each boundary component $\Gamma_{m+k}$ for $k=1,2,\ldots,\ell$. Thus, we can always choose a branch cut of the logarithm function such that the functions $\gamma_{k}(t)$ given by the first formula in (3.10) are periodic for $t\in J_{j}$ for each $j=m+1,m+2,\ldots,m+\ell$.

b) $\ell^{\prime}=\ell-1$. For this case, $\Gamma_{m+\ell}$ is the external boundary component of $G$. Hence, none of the auxiliary points $\alpha,\alpha_{1},\ldots,\alpha_{m^{\prime}}$ is interior to any of the curves $\Gamma_{m+1},\ldots,\Gamma_{m+\ell-1}$. However, all the auxiliary points $\alpha,\alpha_{1},\ldots,\alpha_{m^{\prime}}$ are interior to the curve $\Gamma_{m+\ell}$. Thus, the winding number of the function $\frac{z-\alpha_{k}}{z-\alpha}$ is always zero along each boundary component $\Gamma_{m+k}$ for $k=1,2,\ldots,\ell$. Hence, we can choose a branch cut of the logarithm function such that the functions $\gamma_{k}(t)$ given by the second formula in (3.10) are periodic for $t\in J_{j}$ for each $j=m+1,m+2,\ldots,m+\ell$. For this case, we need to prove also that equation (3.11) holds for the functions $\gamma_{k}(t)$ defined by the second formula in (3.10). Since $\Gamma_{m+\ell}$ is the external boundary component of $G$, we have $m^{\prime}=m$, and by (2.9), we have $\sum_{k=1}^{m^{\prime}}a_{k}=0$. Thus,

[TABLE]

and hence (3.11) holds for the functions $\gamma_{k}(t)$ defined by the second formula in (3.10). ∎

Taking into account (3.11), we rewrite the boundary condition (3.8) as

[TABLE]

where the functions $\gamma_{k}$ are defined by (3.10). Since we are interesting in computing only $u={\rm Re}\,F$, we can assume that $g(\infty)=c$ is real for unbounded $G$ and $g(\alpha)=c$ is real for bounded $G$. We introduce an auxiliary function $f$ defined in $G$ by

[TABLE]

Then $f$ is a single-valued analytic function in $G$ with $f(\infty)=0$ for unbounded $G$. Let $A(t)$ be the complex-valued function defined by [Nas-ETNA]

[TABLE]

Hence the boundary condition (3.12) implies that the function $f$ is a solution of the following Riemann-Hilbert problem

[TABLE]

Observe that solving the Riemann-Hilbert problem (3.15) requires finding the unknown analytic functions $f$ as well as the unknown real constants $a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}$ in the right-hand side of (3.15).

4. The generalized Neumann kernel

The generalized Neumann kernel $N(s,t)$ is defined for $(s,t)\in J\times J$ by [Weg-Nas]

[TABLE]

Closely related to the kernel $N$ is the following kernel $M(s,t)$ defined for $(s,t)\in J\times J$ by [Weg-Nas]

[TABLE]

The kernel $N$ is continuous and the kernel $M$ has a singularity of cotangent type [Weg-Nas].

Let $H$ denote the space of all real-valued Hölder continuous functions on the boundary $\Gamma$. In this paper, for simplicity, if $\phi$ is a real-valued function defined on the boundary $\Gamma$, then we write $\phi(\eta(t))$ as $\phi(t)$. Further, any piecewise constant function $h\in H$ defined by

[TABLE]

with real constants $h_{k}$ for $k=1,\ldots,m+\ell$ will be denoted by

[TABLE]

The integral operators with the kernels $N(s,t)$ and $M(s,t)$ are defined on $H$ by

[TABLE]

The identity operator on $H$ will be denoted by ${\bf I}$. Then, we have the following theorem from [Nas-JMAA1].

Theorem 4.3**.**

For each $k=1,2,\ldots,m^{\prime}$, let the function $\gamma_{k}$ be given by (3.10). Then, there exists a unique real-valued function $\mu_{k}\in H$ and a unique piecewise constant real-valued function $h_{k}=(h_{1,k},h_{2,k},\ldots,h_{m+\ell,k})$ such that

[TABLE]

are boundary values of an analytic function $f_{k}$ in $G$ with $f(\infty)=0$ for unbounded $G$. The function $\mu_{k}$ is the unique solution of the integral equation

[TABLE]

and the function $h_{k}$ is given by

[TABLE]

The integral equation (4.5) been used for computing the conformal map from bounded and unbounded multiply connected domains onto several canonical slit domains, see e.g., [Nas-Siam1, Nas-JMAA1, Nas-ETNA]. The following lemma is needed to prove Theorems 4.9 and 4.22 below.

Lemma 4.7**.**

If $f$ is an analytic function in $G$ with $f(\infty)=0$ for unbounded $G$ such that its boundary values satisfy the boundary condition

[TABLE]

for a piecewise constant real-valued function $\gamma(t)=(c_{1},c_{2},\ldots,c_{m+\ell})$, then $f$ is the zero function and $c_{1}=c_{2}=\cdots=c_{m+\ell}=0$.

Proof.

The solvability of the Riemann-Hilbert problem (4.8) depends on the winding number of the function $A$. For the function $A$ defined in (3.14), the Riemann-Hilbert problem (4.8) is not necessarily solvable [Nas-ETNA]. However, by Theorem 4.3, a unique piecewise constant real-valued function $h(t)=(h_{1},h_{2},\ldots,h_{m+\ell})$ exists such that the Riemann-Hilbert problem

[TABLE]

is uniquely solvable (see also [Nas-ETNA, Weg-Nas]). By the uniqueness of the piecewise constant function $h$ and since the function $\gamma$ is a piecewise constant function, the function $h$ must be given by $h(t)=-\gamma(t)$ since the problem

[TABLE]

will be solvable and has the zero solution $f(z)=0$. ∎

In the remaining part of this section, we shall use Theorem 4.3 to present a method for computing the real constants $a_{1},\ldots,a_{m}$ and hence computing $\mathrm{cap}(C)$ through (2.12). Recall from (2.1) that either $m^{\prime}=m$ or $m^{\prime}=m-1$. These two cases of $m^{\prime}$ will be considered separately in the following two subsections.

4.1. Case I: $m^{\prime}=m$

This case includes the following two subcases:

1. (1)

Both $G$ and $B$ are unbounded (see Figure 1). For this subcase, we have $m^{\prime}=m\geq 2$, $\ell^{\prime}=\ell\geq 0$ (where $B={\mathbb{C}}$ for $\ell=0$), $A$ is given by the first formula in (3.14), and the functions $\gamma_{k}$ for $k=1,2,\ldots,m$ are given by the first formula in (3.10). 2. (2)

Both $G$ and $B$ are bounded (see Figure 2). For this subcase, we have $m^{\prime}=m\geq 2$, $\ell^{\prime}=\ell-1\geq 0$, $\Gamma_{m+\ell}$ is the external boundary component of $G$, $A$ is given by the second formula in (3.14), and the functions $\gamma_{k}$ for $k=1,2,\ldots,m$ are given by the second formula in (3.10).

For these two subcases, all the simply connected domains $G_{1},\ldots,G_{m}$ are bounded (see Figures 1 and 2). In Figures 1 and 2, and in all figures throughout the paper, the boundaries of the domain $B$ are the “dash-dotted” curves and the boundaries of the plates of the condenser are the “solid” curves.

The following theorem provides us with a method for computing the unknown real constants $a_{1},\ldots,a_{m}$. The theorem will be proved using an approach similar to the approach used in proving Theorems 4.2 and 4.3 in [NLS16],

Theorem 4.9**.**

For each $k=1,2,\ldots,m$, let the function $\gamma_{k}$ be defined by (3.10), $\mu_{k}$ be the unique solution of the integral equation (4.5), and the piecewise constant function $h_{k}=(h_{1,k},h_{2,k},\ldots,h_{m+\ell,k})$ be given by (4.6). Then, the boundary values of the function $f$ in (3.15) are given by

[TABLE]

and the $m+\ell+1$ unknown real constants $a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}$ are the components of the unique solution vector of the linear system

[TABLE]

Proof.

Suppose that $f$ is the analytic function in $G$ with $f(\infty)=0$ for unbounded $G$ and satisfies the boundary condition (3.15). Suppose also that $\hat{f}$ is defined in $G$ by

[TABLE]

where $f_{k}$ are as in Theorem 4.3 and the constants $a_{1},\ldots,a_{m}$ satisfy the condition (2.9). Then $\hat{f}$ is analytic in $G$ with $f(\infty)=0$ for unbounded $G$ and the boundary values of $\hat{f}$ satisfy

[TABLE]

Then the function $\Psi$ defined by $\Psi(z)=\hat{f}(z)-f(z)$ is analytic in $G$ with $\Psi(\infty)=0$ for unbounded $G$. Since $m^{\prime}=m$, it follows from (3.15) and (4.13) that

[TABLE]

The right-hand side is a piecewise constant function, and then Lemma 4.7 implies that $\Psi$ is the zero function and hence $f(z)=\hat{f}(z)$. Thus, (4.10) follows from (4.4) and (4.12). Further, since $\Psi$ is the zero function, the right-hand side of (4.14) is also the zero function and hence

[TABLE]

Since, in view of (3.7), $\cos\theta(t)=1$ for $t\in J_{k}$ for $k=1,2,\ldots,m$ and $\cos\theta(t)=0$ for $t\in J_{k}$ for $k=m+1,m+2,\ldots,m+\ell$, then (4.15) and (2.9) imply that the real constants $a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}$ are the components of a solution vector of the linear system (4.11).

To show that the linear system (4.11) has a unique solution, let $[a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}]^{T}$ be a solution to the homogeneous linear system obtained by assuming that the right-hand side of (4.11) is the zero vector. Then, the homogeneous system implies that

[TABLE]

Assume that the functions $f_{k}$ are as in Theorem 4.3 and $\hat{f}$ is defined by (4.12). Hence, in view of (4.13), the boundary values of the function $\hat{f}$ satisfy

[TABLE]

Then, we define a function $\hat{F}$ in $G$ by

[TABLE]

For unbounded $G$, the function $\hat{F}(z)$ can be written as

[TABLE]

Since $\hat{f}(\infty)=0$ and $\sum_{k=1}^{m}a_{k}=0$, we have $\hat{F}(\infty)=0$. Thus, the function $\hat{F}(z)$ is analytic in $G$ for both cases of bounded and unbounded $G$ but it is not necessarily single valued. In view of (3.14), the boundary values of the function $\hat{F}$ satisfy

[TABLE]

Then by (3.11) and (4.17), we have

[TABLE]

which, in view of (3.7), implies that

[TABLE]

By differentiation both sides of (4.19b) with respect to the parameter $t$, we obtain

[TABLE]

Let the real function $u$ be defined for $z\in G\cup\partial G$ by

[TABLE]

Then $u$ is harmonic in $G$. In view of (3.5), we have

[TABLE]

Thus, by (4.19a), (4.20), and (4.21), the boundary values of $u$ satisfy the mixed-boundary condition

[TABLE]

Since the above mixed boundary value problem has a unique solution, it is clear that the unique solution is the constant function $u(z)=-c$ for all $z\in G\cup\partial G$. Thus the real part of $\hat{F}$ is constant for $z\in G$, and hence, by the Cauchy-Riemann equations, $\hat{F}$ is constant in $G$, say equal to $C$. This implies that $\hat{F}(z)=0$ for all $z\in G$ when $G$ is unbounded since $\hat{F}(\infty)=0$. Then, for all $z\in G$, it follows from (4.18) that

[TABLE]

which implies that that $a_{1}=a_{2}=\cdots=a_{m}=0$ since the functions on the right-hand side are single-valued and the function on the left-hand side is multi-valued. Thus, for bounded $G$, we have $(z-\alpha)\hat{f}(z)=C$ for all $z\in G$. By substituting $z=\alpha$, we find $C=0$ and hence $\hat{F}(z)=0$ for all $z\in G\cup\partial G$. Thus for both cases of bounded and unbounded $G$, we have $F(z)=0$ for all $z\in G\cup\partial G$. Hence, it follows from (4.19) that $c=0$ and $\nu_{1}=\nu_{2}=\cdots=\nu_{\ell}=0$. Thus, the homogeneous linear system has only the trivial solution $[a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}]^{T}={\bf 0}$, and hence the matrix of the linear system (4.11) is non-singular. ∎

4.2. Case II: $m^{\prime}=m-1$

For this case, $G$ is a bounded multiply connected domain of connectivity $m+\ell$ with $m\geq 2$ and $B$ is an unbounded multiply connected domain of connectivity $\ell^{\prime}=\ell\geq 0$ (where $B={\mathbb{C}}$ for $\ell=0$). Here, the simply connected domains $G_{1},\ldots,G_{m-1}$ are bounded, the simply connected domain $G_{m}$ is unbounded, and $\Gamma_{m}$ is the external boundary component of $G$ (see Figure 3). Further, $A$ is given by the second formula in (3.14) and the functions $\gamma_{k}$ for $k=1,2,\ldots,m-1$ are given by the first formula in (3.10). For this case, the values of the unknown real constants $a_{1},\ldots,a_{m-1},c,\nu_{1},\ldots,\nu_{\ell}$ can be computed as in the following theorem. Then $a_{m}$ is computed through (2.10).

Theorem 4.22**.**

For each $k=1,2,\ldots,m-1$, let the function $\gamma_{k}$ be defined by (3.10), let $\mu_{k}$ be the unique solution of the integral equation (4.5), and let the piecewise constant function $h_{k}=(h_{1,k},h_{2,k},\ldots,h_{m+\ell,k})$ be given by (4.6). Then, the boundary values of the function $f$ in (3.15) are given by

[TABLE]

and the $m+\ell$ unknown real constants $a_{1},\ldots,a_{m-1},c,\nu_{1},\ldots,\nu_{\ell}$ are the unique solution of the linear system

[TABLE]

Proof.

The theorem can be proved by the same argument as in the proof of Theorem 4.9. ∎

4.3. Computing the capacity $\mathrm{cap}(C)$ and the potential function $u$

By solving the integral equations (4.5) and then solving the linear system (4.11) (or (4.24)), we obtain the real constants $a_{1},\ldots,a_{m}$. Then, we can compute the capacity $\mathrm{cap}(C)$ from (2.12). We can also compute the boundary values of the auxiliary analytic function $f(z)$ through (4.10) or (4.23). Then the values of $f(z)$ at interior points $z\in G$ can be computed by Cauchy’s integral formula. Since $u(z)={\rm Re}\,[F(z)]$, it follows from (2.6) and (3.13) that the function $u(z)$ is given for $z\in G$ by

[TABLE]

4.4. Outline of the algorithm

The method presented in this section for computing the capacity $\mathrm{cap}(C)$ and the potential function $u$ can be summarized in the following algorithm. Steps 10–12 are needed only if it is required to compute the values of the potential function.

Algorithm 4.26**.**

(Computing the capacity $\mathrm{cap}(C)$ and the potential function $u$).**

Parametrize the boundary components $\Gamma_{j}$ by $\eta_{j}(t)$, $t\in[0,2\pi]$, for $j=1,2,\ldots,m+\ell$, where $\Gamma_{j}$ for $j=1,2,\ldots,m$ are the boundaries of the plates $E_{j}$ of the condenser and $\Gamma_{j}$ for $j=m+1,m+2,\ldots,m$ are the boundary components of the domain $B$.

• 2.

If $G$ is bounded and $B$ is unbounded, then we define $m^{\prime}=m-1$ and $\ell^{\prime}=\ell$. For this case, the plates $E_{1},\ldots,E_{m-1}$ are bounded, the plate $E_{m}$ is unbounded and $\Gamma_{m}$ is the external boundary component of $G$.

• 3.

If both domains $B$ and $G$ are bounded, then we define $m^{\prime}=m$ and $\ell^{\prime}=\ell-1$. For this case, the plates $E_{1},\ldots,E_{m}$ are bounded and $\Gamma_{m+\ell}$ is the external boundary component of $G$.

• 4.

If both domains $B$ and $G$ are unbounded, then we define $m^{\prime}=m$ and $\ell^{\prime}=\ell$. For this case, the plates $E_{1},\ldots,E_{m}$ are bounded.

• 5.

Define the functions $A$ by (3.14).

• 6.

Define the functions $\gamma_{k}$ for $k=1,2,\ldots,m^{\prime}$ by (3.10).

• 7.

For $k=1,2,\ldots,m^{\prime}$, compute the function $\mu_{k}$ by solving the integral equation (4.5) and compute the function $h_{k}$ through (4.6).

• 8.

Compute the $m+\ell+1$ real constants $a_{1},\ldots,a_{m},c,\nu_{1},\ldots,\nu_{\ell}$ by solving one of the linear system (4.11) or (4.24). For $m^{\prime}=m-1$, $a_{m}$ is computed through (2.10).

• 9.

Compute the capacity $\mathrm{cap}(C)$ from (2.12).

• 10.

Compute the boundary values of the analytic function $f$ through (4.10) or (4.23).

• 11.

Compute the values of $f(z)$ for $z\in G$ by the Cauchy integral formula.

• 12.

Compute the values of the potential function $u$ by (4.25).

5. Numerical implementation of the algorithm

The main steps in the Algorithm 4.26 are steps 7 and 8. In step 8, the size of the linear system is usually quite small and hence we solve it using MATLAB “backslash” operator. For step 7, the $m^{\prime}$ integral equation (4.5) are solved using the MATLAB function fbie from [Nas-ETNA]. In the function fbie, the integral equations (4.5) is discretized by the Nyström method with the trapezoidal rule [Atk97, Tre-Trap]. The size of the obtained linear system is usually large. So, in the function fbie, the linear system is solved iteratively using the MATLAB function $\mathtt{gmres}$. The matrix-vector multiplication in $\mathtt{gmres}$ is computed in a fast and efficiently way using the MATLAB function $\mathtt{zfmm2dpart}$ from the toolbox $\mathtt{FMMLIB2D}$ [Gre-Gim12]. The function fbie computes also the $m^{\prime}$ piecewise constant functions $h_{k}$ in (4.6).

For domains with smooth boundaries, we use the trapezoidal rule with equidistant nodes. We discretize each interval $J_{k}=[0,2\pi]$, for $k=1,2,\ldots,m+\ell$, by $n$ equidistant nodes $s_{1},\ldots,s_{n}$ where

[TABLE]

and $n$ is an even integer. We write ${\bf s}=[s_{1},\ldots,s_{n}]$. Then, we discretize the parameter domain $J$ by the $m+\ell$ copies of ${\bf s}$,

[TABLE]

This leads to the discretizations

[TABLE]

In MATLAB, these discretized functions are stored in the vectors et, etp, A, gamk, respectively. Then the discretizations vectors muk and hk of the functions $\mu_{k}$ and $h_{k}$ in (4.5) and (4.6) are computed by calling

[muk,hk] = fbie(et,etp,A,gamk,n,iprec,restart,tol,maxit).

In the numerical experiments in the next sections, we choose $\mathtt{iprec}=5$ (the tolerance of the FMM is $0.5\times 10^{-15}$), restart=[ ] (GMRES is used without restart), tol=1e-14 (the tolerance of the GMRES method is $10^{-14}$), and maxit=100 (the maximum number of GMRES iterations is $100$). The values $h_{j,k}$ are then computed by taking arithmetic means:

[TABLE]

These values are used to build the linear system (4.11) or (4.24). Thus, the computational cost of the overall method for computing the capacity $\mathrm{cap}(C)$ is $O(m^{\prime}(m+\ell)n\ln n)$ operations for step (7) and $O((m+\ell)^{3})$ operations for step (8).

For fast and accurate computing of the Cauchy integral formula in step (11), we use the MATLAB function fcau from [Nas-ETNA]. The function fcau is based on using the MATLAB function zfmm2dpart in [Gre-Gim12]. Using the function fcau, the Cauchy integral formula can be computed at $p$ interior points in $O(p+(m+\ell)n)$ operations.

For domains with corners (excluding cusps), the trapezoidal rule with equidistant nodes yields only poor convergence and hence the trapezoidal rule with a graded mesh will be used [Kre90]. Equivalently, we can remove the discontinuity of the derivatives of the solution of the integral equation at the corner points by choosing an appropriate one-to-one function $\sigma:J\to J$. Then we parametrize the boundary $\Gamma$ by $\eta(t)=\hat{\eta}(\sigma(t))$ where $\hat{\eta}$ is any parametrization function of the boundary $\Gamma$ (see [Kre90, LSN17] for more details, the above function $\sigma$ is denoted by $\delta$ in [LSN17]).

The proposed method can be implemented in MATLAB as in the following function capgc.m.

⬇

ction [cap , uz] = capgc(et,etp,alphav,deltav,m,mp,ell,alpha,z)

% Compute the capacity of the generalized condensers (B,E,delta)

%

% Input:

% 1,2) et, etp: parametrization of the boundary and its first derivative

% 3) alphav=[alphav(1),...,alphav(mp)]: alphav(j) is an auxiliary point

% interior to the boundary component \Gamma_j

% 4) deltav=[deltav(1),...,deltav(m)]: deltav(j) is the value of the

% potential function u on \Gamma_j

% 5) m: the number of the closed sets E_k

% 6) mp: mp=m-1 if \Gamma_m is the external boundary component of G,

% o.w., mp=m

% 7) ell: the multiplicity of the domain B (B=C for ell=0)

% 8) alpha: for bounded G, alpha is an auxiliary point in G

% for unbounded G, alpha=inf

% 9) z: a row vector of points in G (if it is required to compute u(z))

%

% Output:

% cap (the capacity of the generalized condensers (C,E,delta)).

% uz (the values of the potential function u(z) if z is given).

%

% Computing the constants \h_{j,k} for j=1,2,...,m+ell and k=1,2,...,mp

ellp = ell ; ellp(abs(alpha)<inf & mp==m)=ell-1;

n=length(et)/(m+ell); tht=zeros(size(et)); tht(m*n+1:end)=pi/2;

if mp==m & ellp==ell

A=exp(-i.*tht);

else

A=exp(-i.*tht).*(et-alpha);

end

for k=1:mp

for j=1:m+ell

    jv = 1+(j-1)*n:j*n;

    if (ellp==ell)

        gamk{k}(jv,1)=real(exp(-i.*tht(jv)).*clog(et(jv)-alphav(k)));

    else

        gamk{k}(jv,1)=real(exp(-i.*tht(jv)).*...

            clog((et(jv)-alphav(k))./(et(jv)-alpha)));

    end

end

[mu{k},h{k}]=fbie(et,etp,A,gamk{k},n,5,[],1e-14,100);

for j=1:m+ell

    hjk(j,k)=mean(h{k}(1+(j-1)*n:j*n));

end

end

% Computing the constants a_k for k=1,2,...,m

mat=hjk; mat(1:m,mp+1)=1; mat(m+1:m+ell,mp+1)=0;

mat(1:m,mp+2:mp+ell+1)=0; mat(m+1:m+ell,mp+2:mp+ell+1)=-eye(ell);

rhs(1:m,1)=deltav; rhs(m+1:m+ell,1)=0;

if mp==m

mat(m+ell+1,1:m)=1; mat(m+ell+1,m+1:m+ell+1)=0; rhs(m+ell+1,1)=0;

end

x=mat\rhs; a=x(1:mp,1); c=x(mp+1);

if mp==m-1

a(m,1)=-sum(a);

end

% Computing the capacity

cap = (2*pi)*sum(deltav(:).*a(:));

% compute the values of the potential function u(z) if z is given

if nargin==9

fet = zeros(size(et)); uz=zeros(size(z));

for k=1:mp

    fet = fet+a(k).*(gamk{k}+h{k}+i.*mu{k})./A;

    uz=uz-a(k)*log(abs(z-alphav(k)));

end

if abs(alpha)<inf

    fz=fcau(et,etp,fet,z);

    uz=uz+c+real((z-alpha).*fz);

else

    fz=fcau(et,etp,fet,z,n,0);

    uz=uz+c+real(fz);

end

end

end

In this paper, computations were performed in MATLAB R2017a on an ASUS Laptop with Intel(R) Core(TM) i7-8750H CPU @2.20GHz, 2208 Mhz, 6 Core(s), 12 Logical Processor(s), and 16GB RAM. The computation times presented in this paper were measured with the MATLAB tic toc commands. All the computer codes of our computations are available in the internet link https://github.com/mmsnasser/gc.

6. Numerical Examples - Regular Condensers

In this section, we shall consider several numerical examples of regular condensers. Some of these examples either have know capacity or have been considered in the literature. So, we can compare the obtained results with the exact capacity or with known capacity computed by other researchers. For such case, we have $\ell=0$ and $\{\delta_{k}\}_{k=1}^{m}$ containing exactly two different numbers which are $1$ and [math].

6.1*.*

Two circles.

In this example, we consider the generalized condenser $C=(B,E,\delta)$ with $B={\mathbb{C}}$ (and hence $\ell=0$), $E=\{E_{1},E_{2}\}$ (and hence $m=2$), and $\delta=\{0,1\}$. The plates of the condenser are given by $E_{k}=\overline{G_{k}}$, $k=1,2$, where $G_{1}=\{z:|z|<1\}$ and $G_{2}=\{z:|z-a|<r\}$ for $r>0$ and a real number $a$ with $a>1+r$. So, for this example, the generalized condenser reduces to a regular condenser, $\ell^{\prime}=\ell=0$, and $m^{\prime}=m=2$. Thus, the field of the condenser, $G$ is the doubly connected domain in the exterior of the two circles $\Gamma_{1}=\{z\,:\,|z|=1\}$ and $\Gamma_{2}=\{z\,:\,|z-a|=r\}$ (see Figure 4 (left) for $a=2$ and $r=0.5$). The exact value of conformal capacity is given by $\mathrm{cap}(G)=2\pi/\log(1/q)$ where $q$ is obtained by solving the following equation [vuo]

[TABLE]

We use the method presented in Section 5 with $n=2^{10}$ to compute approximate values for the capacity for $a=2$ and for several values of $r$ between $0.01$ and $0.99$. The relative errors for the computed values for this case are presented in Figure 4(right). The level curves of the function $u$ for $a=2$ and $r=0.5$ are shown in Figure 4 (left).