# Computing the logarithmic capacity of compact sets having (infinitely) many components with the Charge Simulation Method

We apply the Charge Simulation Method (CSM) in order to compute the logarithmic capacity of compact sets consisting of (infinitely) many "small" components. This application allows to use just a single charge point for each component. The resulting method therefore is significantly more efficient than methods based on discretizations of the boundaries (for example, our own method presented in [Liesen, Sète, Nasser, 2017]), while maintaining a very high level of accuracy. We study properties of the linear algebraic systems that arise in the CSM, and show how these systems can be solved efficiently using preconditioned iterative methods, where the matrix-vector products are computed using the Fast Multipole Method. We illustrate the use of the method on generalized Cantor sets and the Cantor dust.

## Authors

• 5 publications
• 4 publications
• 2 publications
03/03/2018

### A Generalization of the Exponential-Logarithmic Distribution for Reliability and Life Data Analysis

In this paper, we introduce a new two-parameter lifetime distribution, c...
06/16/2020

### Logarithmic Voronoi cells

We study Voronoi cells in the statistical setting by considering preimag...
05/25/2018

Stake systems which issue stakes as well as coins are proposed. Two suba...
02/11/2020

### Smooth Points on Semi-algebraic Sets

Many algorithms for determining properties of real algebraic or semi-alg...
12/15/2021

### A recursive eigenspace computation for the Canonical Polyadic decomposition

The canonical polyadic decomposition (CPD) is a compact decomposition wh...
02/25/2022

### Conformal capacity and polycircular domains

We study numerical conformal mapping of multiply connected planar domain...
06/01/2021

### Quantifying the Similarity of Planetary System Architectures

The planetary systems detected so far already exhibit a wide diversity o...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

As pointed out by Ransford and Rostand, the computation of the logarithmic capacity of compact subsets of the complex plane is in general a “notoriously hard” problem [24, p. 1499]

. This is particularly true for non-connected sets consisting of many, or even infinitely many components. An important and frequently studied example in this context is given by the classical Cantor middle third set and its generalizations, for which no analytic formula for their logarithmic capacity is known. Among the different approaches to computing the logarithmic capacity of these fractal sets are the method of Ransford and Rostand which uses linear programming

[24], a method of Banjai, Embree and Trefethen based on Schwarz-Christoffel mappings (also see [24]), the study of the spectral theory of orthogonal polynomials by Krüger and Simon [13], and our algorithm based on conformal maps onto lemniscatic domains [14]. Further numerical methods for computing the logarithmic capacity of compact sets can be found, for example, in [4, 5, 23, 25].

In this paper we derive and study an alternative method for numerically approximating the logarithmic capacity of compact sets consisting of very many “small” components, which uses the Charge Simulation Method (CSM) [1, 2, 20], also known as the Method of Fundamental Solutions [3, 30]. In the CSM, the Green function of the complement of the given compact set is approximated by a linear combination of logarithmic potentials which depend on so-called charge points inside each component of the set. Solving the (dense) linear algebraic system for the coefficients of the linear combination then yields the desired approximation of the logarithmic capacity.

For the types of sets we consider, choosing one charge point in each component is sufficient. A similar approach has been used in [12]

, where the CSM was used to solve the harmonic image inpainting problem. The linear algebraic system that needs to be solved therefore is significantly smaller than in methods that are based on discretizing the boundaries of each component (for example, our own Boundary Integral Equation (BIE) method presented in

[14]). We solve the linear algebraic systems arising in the CSM iteratively with the GMRES method [26]. In order to speed up the iteration we use the centrosymmetric structure of the system matrices, the Fast Multipole Method [9] for computing matrix-vector products, and a problem-adapted preconditioner. We apply the new method to generalized Cantor sets and the Cantor dust. These sets consist of infinitely many components, and we compute approximations of their logarithmic capacities by extrapolating from computed logarithmic capacities of finite approximations that consist of many components. Timing comparisons with our method from [14] show the computational efficiency, and comparisons of the computed capacity values with other published results demonstrate the high accuracy of the new method.

The paper is organized as follows. In Section 2 we present the necessary background on the logarithmic capacity and the CSM. In Section 3 we study the approximation of the logarithmic capacity for generalized Cantor sets, and in Section 4 we do the same for the Cantor dust. The paper ends with concluding remarks in Section 5.

## 2 Logarithmic capacity and the CSM

Let be a compact set and denote . Suppose that the unbounded connected component of is regular in the sense that it possesses a Green function . Then the logarithmic capacity of , denoted by , satisfies

 cap(E)=limz→∞exp(log|z|−g(z)); (1)

see, e.g., [22, Thm. 5.2.1], [14, Eq. (2.2)], or Szegő’s original article [28]. Recall that the Green function with pole at infinity of is the real-valued function such that

1. is harmonic in and is bounded in a neighborhood of infinity,

2. is continuous in and vanishes on .

We assume that the boundary of the unbounded domain with consists of a finite number of Jordan curves, in which case the Green function with pole at infinity of exists [7, p. 41].

In the CSM, a harmonic function is approximated by a linear combination of fundamental solutions of the Laplace equation, which are logarithmic potentials in the two-dimensional case [3, 30]. In our setting, the Green function is approximated by a function of the form

 h(z)=c+N∑j=1pjlog|z−wj|,z∈G, (2)

where are undetermined real constants, and where are pairwise distinct. The points are called the charge points, and the coefficients the charges [1]. The function in (2) is harmonic in . Since behaves as for , we require

 N∑j=1pj=1.

The coefficients in (2) are usually determined from the boundary condition on , by imposing the condition in a finite set of collocation points on ; see, e.g., [1, 2, 3, 20, 30]. The number of collocation points is usually at least , so that this procedure leads to a square or overdetermined linear algebraic system [30, p. 12]. Below we will use a slightly different approach to obtain a linear algebraic system for the coefficients.

In view of (1), the approximation of yields

 cap(E)=limz→∞exp(log|z|−g(z))≈limz→∞exp(log|z|−h(z))=e−c. (3)

The following result gives an error bound for this approximation.

###### Lemma 2.1.

Let be compact such that the unbounded connected component of is bounded by Jordan curves, and let be the Green function with pole at infinity of . Let be as in (2), then

 |cap(E)−e−c|≤e−c(M+12M2eM), (4)

where and .

###### Proof.

As above, let be the first charge point. The auxiliary functions

 u(z)=g(z)−log|z−w1|,v(z)=h(z)−log|z−w1|,

are continuous in and harmonic in , including at infinity with and . Since , the Möbius transformation maps onto the bounded domain , whose boundary consists of Jordan curves. The functions and are harmonic in and continuous in . Then, by the maximum principle for harmonic functions on bounded domains (see, e.g., [22, Thm. 1.1.8]),

 |(u∘φ−1)(w)−(v∘φ−1)(w)|≤maxω∈∂R|(u∘φ−1)(ω)−(v∘φ−1)(ω)|,w∈¯¯¯¯R.

This yields

 |u(z)−v(z)|≤maxζ∈∂G|u(ζ)−v(ζ)|=maxζ∈∂G|g(ζ)−h(ζ)|=maxζ∈∂G|h(ζ)|,z∈¯¯¯¯G.

The last equality holds since vanishes on . Taking the limit we obtain

. The estimate (

4) follows from and Taylor’s formula. ∎

On the right hand side of (4) we can replace the (in general unknown) quantity by the upper bound given by the maximum of the (known) function on . In this way we obtain a computable upper bound on the error. We will give an example in Section 3.2.

###### Remark 2.2.

We give another interpretation of the approximation in the case . Then is equivalent to . Since the real numbers are not necessarily rational, the set is the exterior of a generalized lemniscate. This is an unbounded domain (with ), and is the Green function with pole at infinity of . In particular, is the logarithmic capacity of the compact set . If has connectivity , it is a lemniscatic domain; see [27, 18] and references therein as well as Walsh’s original article [29]. Otherwise, the connectivity of is lower, and is the exterior of a level curve of the Green function of a lemniscatic domain. Thus, in the method we propose here, the original domain is approximated by the domain , and the method returns as approximation of .

If not all charges are positive, then the set is a canonical domain of the more general form in [29, Thm. 3].

In the following two sections we will apply the CSM in order to compute (approximations of) the logarithmic capacity of generalized Cantor sets (Section 3) and the Cantor dust (Section 4).

## 3 Generalized Cantor sets

In this section we start with setting up the CSM for the generalized Cantor sets. Next we give an analytic example (illustrating this approach and Lemma 2.1), and study the structure and properties of the matrices arising in the CSM. We then show how the resulting linear algebraic systems can be solved iteratively, and finally we present the results of numerical computations of the logarithmic capacity of generalized Cantor sets.

### 3.1 Setting up the CSM

Fix some , let , and define recursively

 Ek:=qEk−1∪(qEk−1+(1−q)),k≥1. (5)

Thus, the set is obtained by removing the middle from each interval of the set ; see Figure 1 for (left) and (right) corresponding to . Then the generalized Cantor set is defined as

 E:=∞⋂k=0Ek. (6)

For we obtain the classical Cantor middle third set. The limiting cases are for and for .

The set consists of intervals , , numbered from left to right. The intervals have same length

 |Ik,1|=|Ik,2|=…=|Ik,m|=qk. (7)

Denote the midpoint of by . Let

 rk=12qk=wk,1, (8)

then

 0

Instead of the sets consisting of closed intervals we will use sets consisting of closed disks in the CSM. More precisely, let

 D0\coloneq{z∈C:∣∣∣z−12∣∣∣≤12},

and define recursively

 Dk:=qDk−1∪(qDk−1+1−q),k≥1.

Then consists of disjoint disks with the centers , , i.e.,

 Dk=m⋃j=1Dk,j;

see Figure 1 for (left) and (right) corresponding to . All of the disks have the same radius from (8). Since , we obtain for all by induction. Moreover, we have the following result.

###### Theorem 3.1.

In the notation established above,

 ∞⋂k=0Dk=Eandcap(E)=limk→∞cap(Dk).
###### Proof.

Clearly, . Next we show that if , then . If , then there exists such that , hence . If , i.e., , then there exists such that , which shows that . Together we obtain . Finally, by [22, Thm. 5.1.3], since are compact and . ∎

Our overall strategy for computing an approximation of , where the set consists of infinitely many components, is to first approximate for reasonably many and large values of  with the CSM. Then we extrapolate from the computed approximations of to obtain an approximation of .

The set is an unbounded multiply connected domain of connectivity with

 ∂Gk=Ck,1∪⋯∪Ck,m,

where is the circle with center and radius for .

As described in Section 2, we approximate , the Green function with pole at infinity of , with the CSM by

 hk(z)=ck+m∑ℓ=1pk,ℓlog|z−wk,ℓ|,z∈Gk, (10)

with

 m∑ℓ=1pk,ℓ=1. (11)

In the CSM for unbounded multiply connected domains, we usually choose many charge points inside each boundary component ,  [2, 20]. However, because is very small for large , we choose only one point inside , which is its center .

###### Remark 3.2.

The fact that we use just a single charge point for each boundary component is an essential difference to other discretization-based methods for computing the logarithmic capacity of sets consisting of many components. For example, our own BIE method presented in [14] “opens up” the intervals of to obtain a compact set of the same capacity, but bounded by smooth Jordan curves. The computation of the logarithmic capacity then is based on discretizing the boundary curves using points on each of them. Consequently, for the same the linear algebraic systems in the method presented here are smaller by a factor of  compared to those in [14]. This is one of the main reasons for the very significant computational savings that we obtain with the new method; see Section 3.5. Using only a single charge point for each boundary component can be justified by the fact that the boundary components become very small when increases. This is illustrated in [12] where, as mentioned in the Introduction, a similar approach has been used. (In terms of the current paper, the derivation of the system matrix in [12] assumes that the value of in the equation (14) is negligible in comparison with for .)

Since vanishes on the boundary , we also require, using the parametrization

 ηk,j(t)=wk,j+rkeit,0≤t≤2π,j=1,2,…,m, (12)

of the circle , that

 0=hk(ηk,j(t)) =ck+m∑ℓ=1pk,ℓlog|wk,j+rkeit−wk,ℓ|, (13) =ck+pk,jlogrk+m∑ℓ=1ℓ≠jpk,ℓlog|wk,j+rkeit−wk,ℓ|,j=1,…,m, (14)

for all . For , the function is harmonic in the disk , hence

 12π∫2π0log|wk,j+rkeit−wk,ℓ|dt=log|wk,j−wk,ℓ|

by the mean value property of harmonic functions; see, e.g., [31, Thm. 4.6.7] or [22, Thm. 1.1.6]. By integrating (14) with respect to the parameter , we thus obtain a linear algebraic system of the form

 ck+pk,jlogrk+m∑ℓ=1ℓ≠jpk,ℓlog|wk,j−wk,ℓ|=0,j=1,2,…,m. (15)

We will now consider (15) for a fixed size , and therefore we will drop the index  in the following for simplicity. We write (15) in the form

 Ap=ce, (16)

where , , and

 A=−⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣logrlog|w1−w2|⋯log|w1−wm−1|log|w1−wm|log|w2−w1|logr⋯log|w2−wm−1|log|w2−wm|⋮⋮⋱⋮⋮log|wm−1−w1|log|wm−1−w2|⋯logrlog|wm−1−wm|log|wm−w1|log|wm−w2|⋯log|wm−wm−1|logr⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (17)

Note that and that is even.

###### Theorem 3.3.

For , the entries of satisfy

 log11−2r≤aij≤log1r,1≤i,j≤m, (18)

and they decay away from the diagonal in each row and column, i.e.,

 ai,1<…ai,i+1>…>ai,m, (19) a1,j<…aj+1,j>…>am,j, (20)

for .

###### Proof.

Since , we have that

 |wi−w1|>|wi−w2|>…>|wi−wi−1|>r,r<|wi−wi+1|<…<|wi−wm|, (21)

which is equivalent to (19). Since is symmetric, (19) is equivalent to (20).

Let . The upper bound in (18) follows from . For the lower bound, notice that by (9), and that if and only if . We then obtain for , which also holds for , since is the largest element in the row. ∎

We will continue under the assumption that is nonsingular, which will be justified below; see Section 3.3 as well as Figure 3 and the corresponding discussion. Note that (11) can be written as , and therefore multiplying (16) from the left with yields or, equivalently,

 c=1eTA−1e. (22)

Thus, in order to compute , which is our approximation of (see (3)), we need to compute (or accurately estimate) the quantity , preferably without explicitly inverting the full matrix . One option is to numerically solve the linear algebraic system

 Ax=e, (23)

and then to compute .

### 3.2 An analytic example: Two disks

In this section, we study the accuracy of the CSM approximation and the tightness of the bounds in Lemma 2.1 on a simple example where an exact formula for the logarithmic capacity is known. We consider the set consisting of the union of the two disks with the centers

 w1=16,w2=56,

and the radius , where . Denote by

 G={z:|z−1/6|>r and |z−5/6|>r}

the complement of in the extended complex plane. Then maps onto the domain , and it follows from [22, Theorem 5.2.3] and [27, Theorem 4.2] (see also [14, Example 4.6]) that

 cap(D1)=cap(ˆGc)=2Kπ√(1/3)2−r2√2L(1+L2)=2K3π√1−9r2√2L(1+L2), (24)

where

 K=K(L2)=∫101√(1−t2)(1−L4t2)dt,L=2ρ∞∏k=1(1+ρ8k1+ρ8k−4)2,

and

 ρ=√1/3+r−√1/3−r√1/3+r+√1/3−r=√1+3r−√1−3r√1+3r+√1−3r=3r1+√1−9r2.

Thus, for we have

 ρ≈3r/2,L≈2ρ≈3r≪1,

hence

 K=K(L2)≈K(0)=π2,

and now (24) implies that

 cap(D1)≈13×1×√2L≈13√6r=√2r3(for r≪1/6). (25)

For the CSM we set up the linear algebraic system (23) with the system matrix

 A=−[logrlog|w1−w2|log|w2−w1|logr]=⎡⎣−logr−log23−log23−logr⎤⎦

from (17), so that

 A−1=1(logr)2−(log23)2⎡⎣−logrlog23log23−logr⎤⎦,

and hence (see (22))

 c=−(logr)2−(log23)22(logr−log23)=−12(logr+log23)=−12log2r3. (26)

Consequently, the CSM estimate is

 cap(D1)≈exp(−c)=exp(12log2r3)=√2r3. (27)

Comparing (25) and (27) shows that the CSM gives very accurate results for . This conclusion is illustrated in Figure 2, where the (blue) solid line shows that the error between the exact capacity (24) and the CSM estimate (27) is for , but increases to for .

We next consider the bound (4) from Lemma 2.1, i.e., the inequality

 |cap(D1)−e−c|≤e−c(M+12M2eM), (28)

where , and

 h(z)=c+p1log|z−w1|+p2log|z−w2|,z∈G.

Here, is given by (26), is given by (24), , and is the solution of the linear algebraic system (16), i.e.,

 ⎡⎣−logr−log23−log23−logr⎤⎦[p1p2]=[cc].

Thus, in view of (26), we have , and hence

 h(z)=12log32r+12log∣∣∣z−16∣∣∣+12log∣∣∣z−56∣∣∣,

which yields

 M≤maxζ∈∂D1|h(ζ)|=12log32+12log(r+23)=12log(32r+1).

Using this upper bound on in (28) we obtain

 |cap(D1)−e−c|≤√r6log(32r+1)(1+14log(32r+1)√32r+1)=:^M. (29)

The values of the computable upper bound on the absolute error are shown by the dotted line in Figure 2. We observe that for larger values of the bound becomes quite tight. Figure 2 also shows the error bound , which we can compute in this example since is known analytically. Clearly, the bound is very tight in this example.

### 3.3 Structure and properties of the system matrices

In Section 3.2 we have considered the case of just two disks, and we were able to invert the -matrix of the linear algebraic system (23) explicitly. For approximating the logarithmic capacity of the generalized Cantor sets by extrapolating from the values for reasonably large , we will have to solve much larger linear algebraic systems, and thus we need to apply more sophisticated techniques for the numerical solution of (23).

Our first observation in this direction is that the matrix in (17) is symmetric as well as centrosymmetric, which means that

 JmAJm=A,whereJm=⎡⎢⎣1\udots1⎤⎥⎦∈Rm×m.

Because of the centrosymmetry, the matrix  can be block-diagonalized with an orthogonal transformation at no additional cost; see [6]. If we partition

 A=[A11A12A21A22],whereA11,A22∈Rm2×m2,

then

and

 B=A11+A12Jm/2,C=A11−A12Jm/2. (30)

Since is symmetric, the matrices are also symmetric. If we partition

 x=[x1x2]

and use the orthogonal decomposition of , then (23) can be transformed via a left-multiplication with into the equivalent system

 [B00C][x1+Jm/2x2x1−Jm/2x2]=[2e0],

where . Since we assume that and thus is nonsingular, the second block row implies that , and hence it remains to solve the system

 By=e, (31)

where , to obtain . In finite precision, we obtain a computed approximate solution of (31), which leads to

 c≈12eT˜y.

Figure 3 shows the -norm condition numbers computed in MATLAB111All computations in this paper are performed in MATLAB R2017a on an ASUS Laptop with Intel Core i7-8750H CPU @ 2.20GHz, 2208 Mhz, 6 Cores, 12 Logical Processors and 16 GB RAM. of and for (i.e., the classical Cantor middle third set) as functions of . We observe that the condition numbers grow linearly in . A similar behavior of the condition numbers can be seen in the Cantor dust example (see Figure 10), and has been observed in [10, Fig. 6], where the matrix in [10, Eq. (36)] has the same form as our matrices and . According to [10, p. 398], such behavior of the condition numbers in the CSM is expected, since the matrices “resemble a discretization of the kernel of a single-layer potential whose inverse is the Laplacian operator”. Analyses of the invertibility of matrices in other applications of the CSM can be found in [19].

From a numerical point of view, solving (31) is clearly preferable over solving (23), since has only half the size of , while both matrices are dense and have essentially the same condition number (cf. Figure 3). We will solve the linear algebraic system (31) using iterative methods that require one matrix-vector product with  in every step; see Section 3.4. Because of the structure of the entries of , this multiplication can be performed using the Fast Multipole Method (FMM) [9]. Using the definition of in (30), each matrix-vector product of the form requires two applications of the FMM. The following result shows that can be written in a form so that only one application of the FMM is required.

###### Lemma 3.4.

The entries of are given by

 bij={−log(2r√zi),i=j,−log|zi−zj|,i≠j,1≤i,j≤m/2,

where for .

###### Proof.

First note that by definition the entries of are given by

 aij={−logr,i=j,−log|wi−wj|,i≠j,1≤i,j≤m/2,

and the entries of are given by

 ^aij=−log|wi−wm/2+j|,1≤i,j≤m/2.

By construction, the , , are real numbers in the interval , which are symmetric about . Further, we have

 wm/2+j=1−q+wj,wm/2+1−j+wj=q,j=1,…,m/2.

Thus,

 ^aij=−log|wi−wj−1+q|,1≤i,j≤m/2,

and hence the entries of the matrix are given by

 ~aij=^ai,m/2+1−j=−log|wi−wm/2+1−j−1+q|=−log|wi+wj−1|,1≤i,j≤m/2.

Finally, by (30), the entries of are given for by

 bii=aij+~aij=−logr−log|2wi−1|=−log|2r(wi−1/2)|,

and for by

 bij=aij+~aij=−log|wi−wj|−log|wi+wj−1|,=−log|(wi−1/2)2−(wj−1/2)2|,

which completes the proof. ∎

From this lemma we have

 B=−⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣log(2r√z1)log|z1−z2|⋯log|z1−zm/2−1|log|z1−zm/2|log|z2−z1|log(2r√z2)⋯log|z2−zm/2−1|log|z2−zm/2|⋮⋮⋱⋮⋮log|zm/2−1−z1|log|zm/2−1−z2|⋯log(2r√zm/2−1)log|zm/2−1−zm/2|log|zm/2−z1|log|zm/2−z2|⋯log|zm/2−zm/2−1|log(2r√zm/2)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦, (32)

and we can easily see that is symmetric (as already mentioned above), but not centrosymmetric.

###### Theorem 3.5.

For , the entries of satisfy

 log1q(1−q)

and they decay away from the diagonal in each row and column, i.e.,

 bi,1<…bi,i+1>…>bi,m/2, (34) b1,j<…bj+1,j>…>bm/2,j, (35)

for .

###### Proof.

In (33), we need , which is equivalent to and thus to .

Decay. Since

 0

we have

 1/4>z1>z2>…>zm/2>0. (37)

This implies

 |zi−z1|>|zi−z