Mobile disks in hyperbolic space and minimization of conformal capacity

TL;DR

This paper investigates the conformal capacity of disk arrangements in hyperbolic space, using numerical methods to find lower bounds and extremal configurations, revealing patterns similar to animal colony heat minimization.

Contribution

It introduces numerical approaches to estimate conformal capacity for complex hyperbolic disk constellations and identifies extremal configurations through capacity minimization.

Findings

Numerical methods effectively estimate capacity bounds.

Extremal configurations feature disks grouped together.

Patterns resemble animal colonies minimizing heat flow.

Abstract

Our focus is to study constellations of disjoint disks in the hyperbolic space, the unit disk equipped with the hyperbolic metric. Each constellation corresponds to a set $E$ which is the union of $m>2$ disks with hyperbolic radii $r_j>0, j=1,...,m$. The centers of the disks are not fixed and hence individual disks of the constellation are allowed to move under the constraints that they do not overlap and their hyperbolic radii remain invariant. Our main objective is to find computational lower bounds for the conformal capacity of a given constellation. The capacity depends on the centers and radii in a very complicated way even in the simplest cases when $m=3$ or $m=4$. In the absence of analytic methods our work is based on numerical simulations using two different numerical methods, the boundary integral equation method and the $hp$-FEM method, resp. Our simulations combine capacity computation with minimization methods and produce extremal cases where the disks of the constellation are grouped next to each other. This resembles the behavior of animal colonies minimizing heat flow in arctic areas.

Tables

Table 1. Table 1. Disks with collinear centers: m 𝑚 m hyperbolic disks E _ ​ j subscript 𝐸 _ 𝑗 E_{\_}j with radii r _ ​ j subscript 𝑟 _ 𝑗 r_{\_}j and centers on the diameter ( − 1 , 1 ) 1 1 (-1,1) with ∑ _ j = 1 m ​ 2 ​ r _ ​ j = d _ ​ 1 = ρ _ ​ 𝔹 2 ​ ( − 0.6 , 0.6 ) subscript _ 𝑗 superscript 1 𝑚 2 subscript 𝑟 _ 𝑗 subscript 𝑑 _ 1 subscript 𝜌 _ superscript 𝔹 2 0.6 0.6 \sum_{\_}{j=1}^{m}2r_{\_}j=d_{\_}1=\rho_{\_}{\mathbb{B}^{2}}(-0.6,0.6) . Conjectured lower and upper bounds of the capacity cap ​ ( 𝔹 , ∪ E _ ​ j ) cap 𝔹 subscript 𝐸 _ 𝑗 {\rm cap}(\mathbb{B},\cup E_{\_}j) .

$m$	Lower	Upper
$2$	$8.515312094751020$	$11.463763614692954$
$3$	$7.450131756754710$	$12.744594178229441$
$4$	$7.017838565418236$	$14.282099489357595$

Table 2. Table 2. Computed values of cap ​ ( 𝔹 2 , E ) cap superscript 𝔹 2 𝐸 {\rm cap}(\mathbb{B}^{2},E) when m = 4 𝑚 4 m=4 for a constellation with disk radii (from left to right) r _ ​ 1 = 0.15 ​ d _ ​ 1 / 2 subscript 𝑟 _ 1 0.15 subscript 𝑑 _ 1 2 r_{\_}1=0.15d_{\_}1/2 , r _ ​ 2 = 0.35 ​ d _ ​ 1 / 2 subscript 𝑟 _ 2 0.35 subscript 𝑑 _ 1 2 r_{\_}2=0.35d_{\_}1/2 , r _ ​ 3 = 0.20 ​ d _ ​ 1 / 2 subscript 𝑟 _ 3 0.20 subscript 𝑑 _ 1 2 r_{\_}3=0.20d_{\_}1/2 , and r _ ​ 4 = d _ ​ 1 / 2 − ( r _ ​ 1 + r _ ​ 2 + r _ ​ 3 ) subscript 𝑟 _ 4 subscript 𝑑 _ 1 2 subscript 𝑟 _ 1 subscript 𝑟 _ 2 subscript 𝑟 _ 3 r_{\_}4=d_{\_}1/2-(r_{\_}1+r_{\_}2+r_{\_}3) where d _ ​ 1 = ρ _ ​ 𝔹 2 ​ ( − 0.6 , 0.6 ) subscript 𝑑 _ 1 subscript 𝜌 _ superscript 𝔹 2 0.6 0.6 d_{\_}1=\rho_{\_}{\mathbb{B}^{2}}(-0.6,0.6) . The centers on the diameter ( − 1 , 1 ) 1 1 (-1,1) as a function of the hyperbolic distance d 𝑑 d between disks, i.e., c _ ​ 1 = − th ⁡ ( ( r _ ​ 2 + r _ ​ 1 + d ) / 2 ) subscript 𝑐 _ 1 th subscript 𝑟 _ 2 subscript 𝑟 _ 1 𝑑 2 c_{\_}1=-\operatorname{\mathrm{th}}((r_{\_}2+r_{\_}1+d)/2) , c _ ​ 2 = 0 subscript 𝑐 _ 2 0 c_{\_}2=0 , c _ ​ 3 = th ⁡ ( ( r _ ​ 2 + r _ ​ 3 + d ) / 2 ) subscript 𝑐 _ 3 th subscript 𝑟 _ 2 subscript 𝑟 _ 3 𝑑 2 c_{\_}3=\operatorname{\mathrm{th}}((r_{\_}2+r_{\_}3+d)/2) , c _ ​ 4 = − th ⁡ ( ( r _ ​ 2 + 2 ​ r _ ​ 3 + r _ ​ 4 + 2 ​ d ) / 2 ) subscript 𝑐 _ 4 th subscript 𝑟 _ 2 2 subscript 𝑟 _ 3 subscript 𝑟 _ 4 2 𝑑 2 c_{\_}4=-\operatorname{\mathrm{th}}((r_{\_}2+2r_{\_}3+r_{\_}4+2d)/2) .

$d$	FEM	BIE	Agreement
$0.00$	$7.017838565413617$	—	—
$0.05$	$7.230698262298420$	$7.230698262298405$	$1.51\times {10}^{-14}$
$0.10$	$7.442082617728579$	$7.442082617728490$	$8.88\times {10}^{-14}$
$0.15$	$7.651760366696882$	$7.651760366696745$	$1.37\times {10}^{-13}$
$0.20$	$7.859490827905997$	$7.859490827905935$	$6.22\times {10}^{-14}$
$0.25$	$8.064996233395842$	$8.064996233395734$	$1.08\times {10}^{-13}$
$0.30$	$8.267972932727597$	$8.267972932727497$	$9.95\times {10}^{-14}$

Table 3. Table 3. Permutations of contiguous constellations. E _ ​ D subscript 𝐸 _ 𝐷 E_{\_}D with centers on the segment ( − 1 , 1 ) 1 1 (-1,1) and ( A , B , C , D ) = ( 1 / 2 , 2 / 5 , 1 / 4 , 1 / 5 ) 𝐴 𝐵 𝐶 𝐷 1 2 2 5 1 4 1 5 (A,B,C,D)=(1/2,2/5,1/4,1/5) . E _ ​ C subscript 𝐸 _ 𝐶 E_{\_}C with centers on the circle { z : | z | = 1 / 2 } conditional-set 𝑧 𝑧 1 2 \{z:|z|=1/2\} and ( A , B , C , D ) = ( 1 / 2 , 1 / 3 , 1 / 4 , 1 / 5 ) 𝐴 𝐵 𝐶 𝐷 1 2 1 3 1 4 1 5 (A,B,C,D)=(1/2,1/3,1/4,1/5) .

Case	$r_{\_}1$	$r_{\_}2$	$r_{\_}3$	$r_{\_}4$	$\mathrm{cap}C\left(E_{\_}D\right)$	$\mathrm{cap}C\left(E_{\_}C\right)$
1	D	B	A	C	6.781488018927628	6.451424010111881
2	D	A	B	C	6.788910565780309	6.455800945561348
3	D	C	A	B	6.843774515059010	6.475070264106950
4	C	D	A	B	6.882473842468833	6.485425869048534
5	A	B	C	D	6.890544149275032	6.496389476635198
6	B	C	A	D	6.897202225461369	6.500210100051595
7	C	A	D	B	6.919626376828870	6.520197932005349
8	A	B	D	C	6.928074481413122	6.523073055329720
9	A	C	B	D	6.932436180755356	6.542555705939787
10	C	B	D	A	6.962814943144452	6.542981227003898
11	A	C	D	B	7.053764008325471	6.575258877036491
12	A	D	C	B	7.055565195334228	6.576332514877286

Table 4. Table 4. Hyperbolic radii used in Figure  6 .

Case	$r_{\_}1$	$r_{\_}2$	$r_{\_}3$	$r_{\_}4$
1	0.4	0.2	0.5	0.25
2	0.2	0.5	0.3	0.5
3	0.5	0.5	0.5	0.2
4	0.2	0.7	0.4	0.1

Table 5. Table 5. Solution times for the minimization process when one disk is fixed and three disks are mobile. Number of steps is number of iterations in the interior-point algorithm. Number of evaluations is the total number of solves performed during the minimization.

Method	Discretization	Time	Number of steps	Number of evaluations
BIE	$n=2^4$	472.9	151	1455
	$n=2^7$	85.6	24	192
	$n=2^{10}$	150.7	24	192
$hp$-FEM	$p=4$	39600	202	39568
	$p=6$	11100	37	7494
	$p=8$	9100	20	4150

Equations

$$\mathrm{cap}({\mathbb{B}}^2,{\cup }_{\_}{j=1}^m{E}_{\_}j)$$

$${\text{sh}}^2\frac{{\rho }_{\_}{\mathbb{B}}^2(x,y)}{2}$$

$$m_{\_}H(x,y)=\frac{y\left(1-|x{|}^2\right)+x\left(1-|y{|}^2\right)}{1-|x{|}^2|y{|}^2+A[x,y]\sqrt{(1-|x{|}^2)(1-|y{|}^2)}}$$

$$B_{\_}\rho (x,M)=\{z\in {\mathbb{B}}^2:{\rho }_{\_}{\mathbb{B}}^2(x,z)<M\}$$

$$\{\begin{array}{ll}{B}_{\_}\rho (x,M)=B^2(y,r),& \\ \text{vskip}5.0pt& \\ y=\frac{x(1-t^2)}{1-|x{|}^2{t}^2},r=\frac{(1-|x{|}^2)t}{1-|x{|}^2{t}^2},t=\mathrm{th}(M/2),& \end{array}$$

$$B_{\_}\rho (0,M)=B^2(0,\mathrm{th}(M/2)).$$

$$\{\begin{array}{ll}{B}^2(y,r)=B_{\_}\rho (x,M),& \\ \text{vskip}5.0pt& \\ x=ty/|y|,M={\rho }_{\_}{\mathbb{B}}^2(x,z),t=m_{\_}H(|y|-r,|y|+r),& \end{array}$$

$$\mathrm{cap}(G,E)=\underset{\_}{\inf }u\in A{\int }_{\_}G|\nabla u{|}^{2}dm,$$

$$\mathrm{cap}(G,E)=\underset{\_}{\iint }\Omega |\nabla u{|}^{2}dm.$$

$$\mathrm{cap}(G,E)=\mathsf{M}(\Delta (E,\partial G;G)),$$

$$\underset{\_}{\inf }v\Vert u_{\_}0-v{\Vert }_{\_}{H}^1(\Omega )\le C\exp (-b\sqrt[3]{N}),$$

$$a(\hat{u},v)=l(v)(=a(u,v))(\forall v\in V),$$

$$a(\epsilon ,v)=l(v)-a(\hat{u},v)(=a(u-\hat{u},v))(\forall v\in W).$$

$$\frac{c}{C}\Vert \epsilon {\Vert }_{\_}1\le \Vert u-\hat{u}{\Vert }_{\_}1\le K\left(\Vert \epsilon {\Vert }_{\_}1+\mathrm{osc}(R,r,\mathcal{T})\right),$$

$$\eta(t)=\left\{\begin{array}[]{cc}\eta_{\_}0(t),&t\in J_{\_}0,\\ \eta_{\_}1(t),&t\in J_{\_}1,\\ \vdots\\ \eta_{\_}m(t),&t\in J_{\_}m.\\ \end{array}\right.$$

$$A(t)=\eta (t)-\alpha ,$$

$$N(s,t):=\frac{1}{\pi }\mathrm{Im}\left(\frac{A(s)}{A(t)}\frac{{\eta }^{\prime }(t)}{\eta (t)-\eta (s)}\right).$$

$$M(s,t):=\frac{1}{\pi }\mathrm{Re}\left(\frac{A(s)}{A(t)}\frac{{\eta }^{\prime }(t)}{\eta (t)-\eta (s)}\right),(s,t)\in J\times J.$$

$${\gamma }_{\_}k(t)=\log |\eta (t)-z_{\_}k|,$$

$${\mu }_{\_}k-\mathbf{N}{\mu }_{\_}k=-\mathbf{M}{\gamma }_{\_}k,$$

$$h_{\_}k=[\mathbf{M}{\mu }_{\_}k-(\mathbf{I}-\mathbf{N}){\gamma }_{\_}k]/2.$$

$$\mathtt{et}=[{\eta }_{\_}0(\mathbf{s}),{\eta }_{\_}1(\mathbf{s}),\dots ,{\eta }_{\_}m(\mathbf{s}){]}^T,\mathtt{etp}=[{\eta }_{\_}^{\prime }0(\mathbf{s}),{\eta }_{\_}^{\prime }1(\mathbf{s}),\dots ,{\eta }_{\_}^{\prime }m(\mathbf{s}){]}^T.$$

$$[\mathtt{muk},\mathtt{hk}]=\mathtt{fbie}(\mathtt{et},\mathtt{etp},\mathtt{A},\mathtt{gamk},\mathtt{n},5,[],1e-14,100),$$

$$h_{\_}j,k=\frac{1}{n}\sum _{\_}{i=1+jn}^{(j+1)n}{\mathtt{hk}}_{\_}i,j=0,1,\dots ,m,k=1,\dots ,m.$$

$$\left[\begin{array}[]{ccccc}h_{\_}{0,1}&h_{\_}{0,2}&\cdots&h_{\_}{0,m}&1\\ h_{\_}{1,1}&h_{\_}{1,2}&\cdots&h_{\_}{1,m}&1\\ \vdots&\vdots&\ddots&\vdots&\vdots\\ h_{\_}{m,1}&h_{\_}{m,2}&\cdots&h_{\_}{m,m}&1\\ \end{array}\right]\left[\begin{array}[]{c}a_{\_}1\\ a_{\_}2\\ \vdots\\ a_{\_}{m}\\ c\end{array}\right]=\left[\begin{array}[]{c}0\\ 1\\ \vdots\\ 1\end{array}\right].$$

$$\mathrm{cap}(\mathbb{B},E)=2\pi \sum _{\_}{k=1}^m{a}_{\_}k.$$

$$\underset{\_}{\min }E$$

$$E_{\_}i\cap E_{\_}j$$

$$E_{\_}j$$

$$W=B^2(\mathrm{i}\tau ,\sqrt{1+{\tau }^2})\cap B^2(-\mathrm{i}\tau ,\sqrt{1+{\tau }^2}),\tau >0,$$

$$\mathrm{cap}(\mathbb{B},{\cup }_{\_}{j=1}^m{E}_{\_}j)\approx \sum _{\_}{j=1}^m\mathrm{cap}(\mathbb{B},E_{\_}j),$$

Full text

1. Introduction

Many extremal problems of physics, exact sciences, and mathematics have solutions which exhibit varying degree of symmetry. A typical situation is to minimize or maximize a set functional of a planar set under the constraint that some other functional is constant. The classical isoperimetric problem [31] is an example. Here one maximizes the area of a planar set given its perimeter and the extremal domain is the disk. G. Pólya and G. Szegö [31] initiated a systematic study of a large class of isoperimetric type problems of mathematical physics for domain functionals such as moment of inertia, principal frequency, torsional rigidity, and, in particular, capacities of condensers. Certain geometric transformations, known under the general name “symmetrization” have the property that they decrease the value of domain functionals and thus can give hints about the extremal configuration of isoperimetric problems [3, 8]. We study here new types of transformations which decrease the value of conformal capacity.

In a very interesting recent paper, A. Solynin [34] describes capacity problems, motivated by the behavior of herds of arctic animals which keep close together to minimize the total loss of heat of the herd or to defend against predators (see figures in [34]). Such a herd behavior seems to suggest the heuristic idea that “minimization of herd’s outer perimeter” minimizes the loss of heat or danger from predators. This kind of extremal problem can be classified as special type of isoperimetric problem. As an illustration of the connection between the kind of transformations we are interested in and the observed behavior in nature, see Figure 1.

In a recent paper [29], isoperimetric inequalities in hyperbolic geometry were applied to estimate the conformal capacity of condensers of the form $(\mathbb{B}^{2},E)$ where $E$ is a union of finitely many disjoint closed disks $E_{\_}j,j=1,...,m,$ in the unit disk $\mathbb{B}^{2}$. Thus $E$ is a constellation of disks. Gehring’s lower bound [9] (see also [29]) is given by condensers of the form $(\mathbb{B}^{2},E^{*})$ where $E^{*}$ is a disk with the hyperbolic area equal to that of $\cup_{\_}{j=1}^{m}E_{\_}j.$ Further recent investigations of condenser capacity in the framework of hyperbolic geometry include [28, 26, 27], where pointers to earlier work can be found. It should be noticed that due to the conformal invariance of the conformal capacity, the hyperbolic geometry provides the natural setup for this study.

We continue here this work and our goal is to analyse extremal cases of the aforementioned capacity and how the capacity depends on the geometry of the disk constellation. The constraint that the disks do not overlap leads to problems of combinatorial geometry. Some examples of such geometric problems, related to this work and the herd behavior mentioned above, are Descartes’ problem of four circles with each circle tangent to three circles, Apollonian circle packing, and Soddy’s “complex kiss precise” problem for configurations of mutually tangent circles [21]. Combinatorial geometry extremal problems motivated by biochemistry research and drug development are described in [23]. A very interesting discussion of many topics of combinatorial geometry including packing problems is given in the encyclopedic work of M. Berger [6]. The three dimensional case is much more difficult than the planar case and it is the subject of the extensive review paper [20] where topics range from optimal packing of spheres to constrained motion of small spheres on the surface of the unit sphere. For an extensive survey of potential theoretic extremal problems see [7].

Analysing the extremal cases of the lower bound for

[TABLE]

for a constellation of disjoint hyperbolic disks $E_{\_}j$ seems to be very difficult even in the simplest cases $m=3,4.$ Therefore we consider this problem in special cases such as the case when the circle centers are at the same distance from origin or analyse constrained motion of one circle along three other fixed circles (see Figure 1). Simulations indicate that several constellations yield local minima of the capacity. Throughout, the hyperbolic geometry provides the natural geometric framework for our study, because of the conformal invariance of the capacity. We use two numerical methods for computing the capacity, the $hp$-FEM and the boundary integral equation (BIE) method. The numerical results lead to a number of conjectures and improved bounds. Indeed, the existing lower bound for constellations considered here is improved of the order of 10% for disks of unit hyperbolic radius. Moreover, the asymptotic nature of the theoretical lower bound as the hyperbolic radii $r_{\_}j\to\infty$ is easily understood in the context of hyperbolic geometry.

In modern physics, in particular in condensed matter physics, there has been a lot of interest in geometric settings with negative curvature [19, 22], that is, exactly our natural setup. The purpose of this paper is also to show how computations can be formulated and performed in both Euclidean and hyperbolic geometries, even with the possibility of moving from one to another. This is highlighted in the last section where the optimal configurations in hyperbolic geometry are found by successive transformations to a Euclidean coordinate system employed in the optimization routines. For information about potential theory and its applications, see [7, 31, 32, 36].

The contents are organized into sections as follows. Section 2 contains the key facts about hyperbolic geometry, including the transformation formulae from Euclidean disks to Poincaré disks and back. Section 3 covers the preliminary notations of conformal capacity, collected from various sources, e.g. from [4, 8, 10, 11, 17, 16]. These are the cornerstones of the geometric setup of the computations in the sequel. Section 3 also provides an overview of the $hp$-FEM [15, 14] adjusted to the present computational tasks, our second computational work horse, the BIE method [25, 28], and the interior-point method used in optimization. The numerical experiments are discussed in Sections 4 and 5. In Section 4 the selected configurations have been designed a priori, with the goal of forming an understanding of the identifiable geometric features of the minimal capacity configurations. In Section 5 that understanding is challenged by searching for the minimal capacity configurations using numerical optimization starting with random initial configurations. Finally, the conclusions are drawn in Section 6.

2. From Euclidean Disk to Poincare and Back

In this section the central transformation formulae collected from various sources are presented. In Figure 2 different properties of geometry on the Poincaré disk have been illustrated. In particular, the facts that for all $\epsilon>0$, $M>0$ there are hyperbolic disks with radii $M$ but Euclidean diameter $<\epsilon$ and hyperbolic disks with equal radii have different Euclidean radii depending on their location are important for our discussion below.

For a point $x\in\mathbb{R}^{n}$ and a radius $r>0$, define an open Euclidean ball $B^{n}(x,r)=\{y\in\mathbb{R}^{n}\text{ }|\text{ }|x-y|<r\}$ and its boundary sphere $S^{n-1}(x,r)=\{y\in\mathbb{R}^{n}\text{ }|\text{ }|x-y|=r\}$. For the unit ball and sphere, we use the simplified notations $\mathbb{B}^{n}=B^{n}(0,1)$ and $S^{n-1}=S^{n-1}(0,1)$. The segment joining two points $x,y\in\mathbb{R}^{n}$ is denoted $[x,y].$

Define the hyperbolic metric in the Poincaré unit disk $\mathbb{B}^{2}$ as in [4], [5, (2.8) p. 15]

[TABLE]

We use the notation $\operatorname{\mathrm{sh}}$ and $\,\mathrm{arsh}\,$ for the hyperbolic sine and its inverse, respectively, and similarly, $\operatorname{\mathrm{th}}$ and $\,\mathrm{arth}\,$ for the hyperbolic tangent and its inverse. The hyperbolic midpoint of $x,y\in\mathbb{B}^{2}$ is given by [37]

[TABLE]

where $A[x,y]=\sqrt{|x-y|^{2}+(1-|x|^{2})(1-|y|^{2})}$. We use the notation

[TABLE]

for the hyperbolic disk centered at $x\in\mathbb{B}^{2}$ with radius $M>0\,.$ It is a basic fact that they are Euclidean disks with the center and radius given by [16, p.56, (4.20)]

[TABLE]

Note the special case $x=0$,

[TABLE]

Conversely, the Euclidean disks can be considered as hyperbolic ones by [37]

[TABLE]

Lemma 2.5** ([4, Thm 7.2.2, p. 132]).**

The area of a hyperbolic disc of radius $r$ is $4\pi\operatorname{\mathrm{sh}}^{2}(r/2)$ and the length of a hyperbolic circle of radius $r$ is $2\pi\operatorname{\mathrm{sh}}(r)$.

3. Conformal Capacity and Numerical Methods

A condenser is a pair $(G,E)$, where $G\subset\mathbb{B}^{2}$ is a domain and $E$ is a compact non-empty subset of $G$. The * conformal capacity* of this condenser is defined as [8, 10, 11, 16, 17]

[TABLE]

where $A$ is the class of $C^{\infty}_{\_}0(G)$ functions $u:G\to[0,\infty)$ with $u(x)\geq 1$ for all $x\in E$ and $dm$ is the $2$-dimensional Lebesgue measure. In this paper we assume that $G=\mathbb{B}^{2}$ is the unit disk and $E=\cup_{\_}{j=1}^{m}E_{\_}j$ where $E_{\_}1,\ldots,E_{\_}m$ are $m$ closed disjoint disks in the unit disk. Hence $\Omega=G\backslash E$ is a multiply connected circular domain of connectivity $m+1$. In this case, the infimum is attained by a function $u$ which is harmonic in $\Omega$ and satisfies the boundary conditions $u=0$ on $\partial G$ and $u=1$ on $\partial E$ [8]. The capacity can be expressed in terms of this extremal function as

[TABLE]

The conformal capacity of a condenser is one of the key notions of potential theory of elliptic partial differential equations [17, 11] and it has numerous applications to geometric function theory, both in the plane and in higher dimensions, [8, 10, 11, 16, 17]. Numerous variants of the definition (3.1) of capacity are given in [10, 11]. First, the family $A$ may be replaced by several other families by [10, Lemma 5.21, p. 161]. Furthermore,

[TABLE]

where $\Delta(E,\partial G;G)$ is the family of all curves joining $E$ with the boundary $\partial G$ in the domain $G$ and $\mathsf{M}$ stands for the modulus of a curve family [10, Thm 5.23, p. 164]. For the basic facts about capacities and moduli, the reader is referred to [10, 11, 16, 17].

3.1. Numerical Methods

In this section the numerical methods used in the numerical experiments are briefly described. The capacities are computed using the $hp$-version of the finite element method (FEM) and the boundary integral equation with the generalized Neumann kernel method (BIE). The minimization problems are computed using the interior-point method as implemented in MATLAB and Mathematica.

Since the Dirichlet problem (3.1) is one of the primary numerical model problems, any standard solution technique can be viewed as having been validated. Verification of the results is discussed in connection with one of the numerical experiments below.

3.1.1. $hp$-FEM

What is of particular interest in the context of this paper is that the $hp$-FEM allows for large curved elements without significant loss of accuracy. Since the number of elements can be kept relatively low given that additional refinement can always be added via elementwise polynomial degree, variation in the boundary can be addressed directly at the level of the boundary representation in some exact parametric form. This is illustrated in Figure 3.

The following theorem due to Babuška and Guo [2] sets the limit to the rate of convergence. Notice that construction of the appropriate spaces is technical. For rigorous treatment of the theory involved see Schwab [33] and references therein.

Theorem 3.4**.**

Let $\Omega\subset\mathbb{R}^{2}$ be a polygon, $v$ the FEM-solution of (3.1), and let the weak solution $u_{\_}0$ be in a suitable countably normed space where the derivatives of arbitrarily high order are controlled. Then

[TABLE]

where $C$ and $b$ are independent of $N$, the number of degrees of freedom. Here $v$ is computed on a proper geometric mesh, where the order of an individual element is set to be its element graph distance to the nearest singularity. (The result also holds for meshes with constant polynomial degree.)

Consider the abstract problem setting with $u$ defined on the standard piecewise polynomial finite element space on some discretization $\mathcal{T}$ of the computational domain $\Omega$. Assuming that the exact solution $u\in H_{\_}0^{1}(D)$ has finite energy, we arrive at the approximation problem: Find $\hat{u}\in V$ such that

[TABLE]

where $a(\,\cdot\,,\,\cdot\,)$ and $l(\,\cdot\,)$, are the bilinear form and the load potential, respectively. Additional degrees of freedom can be introduced by enriching the space $V$. This is accomplished via introduction of an auxiliary subspace or “error space” $W\subset H_{\_}0^{1}(D)$ such that $V\cap W=\{0\}$. We can then define the error problem: Find $\varepsilon\in W$ such that

[TABLE]

This can be interpreted as a projection of the residual to the auxiliary space.

The main result on this kind of estimators for the Dirichlet problem (3.1) is the following theorem.

Theorem 3.7** ([14]).**

There is a constant $K$ depending only on the dimension $d$, polynomial degree $p$, continuity and coercivity constants $C$ and $c$, and the shape-regularity of the triangulation $\mathcal{T}$ such that

[TABLE]

where the residual oscillation depends on the volumetric and face residuals $R$ and $r$, and the triangulation $\mathcal{T}$.

3.1.2. BIE method

We review a BIE method from [28] for computing the capacity $\operatorname{\mathrm{cap}}(\mathbb{B},E)$. The method is based on the BIE with the generalized Neumann kernel. The domains considered in this paper are circular domains, i.e., domains whose boundary components are circles. The external boundary is the unit circle, denoted by $C_{\_}0$, is parametrized by $\eta_{\_}0(t)=e^{\mathrm{i}t}$ for $t\in J_{\_}0=[0,2\pi]$. The inner circles $C_{\_}j$ are parametrized by $\eta_{\_}j(t)=z_{\_}j+r_{\_}je^{-\mathrm{i}t}$, $t\in J_{\_}j=[0,2\pi]$, for $j=1,2,\ldots,m$, where $z_{\_}j$ is the center of the circle $C_{\_}j$ and $r_{\_}j$ is its radius. Let $J$ be the disjoint union of the $m+1$ intervals $J_{\_}j=[0,2\pi]$, $j=0,1,\ldots,m$. We define a parametrization of the whole boundary $C=\cup_{\_}{j=0}^{m}C_{\_}j$ on $J$ by (see [25] for the details)

[TABLE]

With the parametrization $\eta(t)$ of the whole boundary $C$, we define a complex function $A$ by

[TABLE]

where $\alpha$ is a given point in the domain $G$. The generalized Neumann kernel $N(s,t)$ is defined for $(s,t)\in J\times J$ by

[TABLE]

We define also the following kernel

[TABLE]

The kernel $N(s,t)$ is continuous and the kernel $M(s,t)$ is singular where the singular part involves the cotangent function. Hence, the integral operator ${\bf N}$ with the kernel $N(s,t)$ is compact and the integral operator ${\bf M}$ with the kernel $M(s,t)$ is singular. Further details can be found in [38].

For each $k=1,2,\ldots,m$, let the function $\gamma_{\_}k$ be defined by

[TABLE]

let $\mu_{\_}k$ be the unique solution of the BIE

[TABLE]

and let the piecewise constant function $h_{\_}k=(h_{\_}{0,k},h_{\_}{1,k},\ldots,h_{\_}{m,k})$ be given by

[TABLE]

For each $k=1,2,\ldots,m$, the solution $\mu_{\_}k$ of the BIE (3.12) and the piecewise constant function $h_{\_}k$ in (3.13) will be computed using the MATLAB fbie from [25]. In the function fbie, the integral equation (3.12) is solved using the Nyström method with the trapezoidal rule. Solving the integral equation is then reduced to solving an $(m+1)n\times(m+1)n$ linear system which is solved by the MATLAB function gmres. The matrix-vector product in gmres is computed by the MATLAB function zfmm2dpart from the MATLAB toolbox FMMLIB2D [12]. To use the MATLAB function fbie, we define a vector ${\bf s}=[s_{\_}1,\ldots,s_{\_}n]$ where $s_{\_}k=2(k-1)\pi/n$, $k=1,\ldots,n$, and $n$ is a given even positive integer. Then we compute the $(m+1)n\times 1$ discretization vectors et and etp of the parametrization $\eta(t)$ of the boundary $C$ and its derivative $\eta^{\prime}(t)$ by

[TABLE]

We also discretize the functions $A(t)$ and $\gamma_{\_}k(t)$ by ${\tt A}={\tt et}-\alpha$ and ${\tt gamk}=\gamma_{\_}k({\tt et})$, $k=1,\ldots,m$. Then we compute $(m+1)n\times 1$ approximate discretizations muk and hk of the functions $\mu_{\_}k(t)$ and $h_{\_}k(t)$ by calling

[TABLE]

i.e., the tolerance of the FMM is $0.5\times 10^{-15}$, the GMRES is used without restart, the tolerance of the GMRES method is $10^{-14}$ and the maximal number of GMRES iterations is $100$.

By computing the $(m+1)n\times 1$ vector hk, we obtain approximate discretizations of the piecewise constant function $h_{\_}k=(h_{\_}{0,k},h_{\_}{1,k},\ldots,h_{\_}{m,k})$ in (3.13). Note that, for $k=1,\ldots,m$, the constant $h_{\_}{j,k}$ is the value of the function $h_{\_}k$ on the boundary component $\Gamma_{\_}j$. We approximate the values of the real constants $h_{\_}{j,k}$ by taking arithmetic means

[TABLE]

The values of the $m$ real constants $a_{\_}1,\ldots,a_{\_}{m}$ are then approximated by solving the $(m+1)\times(m+1)$ linear system [28]

[TABLE]

Since $m+1$ is the number of boundary components of the domain $\Omega=G\setminus E$, we can assume that $m$ is small and solve the linear system (3.14) using the Gauss elimination method. By solving the linear system, the capacity $\operatorname{\mathrm{cap}}(\mathbb{B},E)$ will be computed by [28, Eq. (3.9)]

[TABLE]

In this paper, the boundary components of the domain $\Omega$ are circles. Thus, the integrands in (3.12) and (3.13) will be $2\pi$-periodic functions, and can be extended holomorphically to some parallel strip $|{\rm Im}\,t|<\sigma$ in the complex plane. Hence, the trapezoidal rule will then converge exponentially with $O(e^{-\sigma n})$ [35] when it is used to discretize the integrals in (3.12) and (3.13). The numerical solution of the integral equation will converge with a similar rate of convergence [1, p. 322] (see Figure 5 (right) below).

3.1.3. Nonlinear Optimization: Interior-Point Method

The two methods outlined above are combined with a numerical optimization routine in the last set of numerical experiments below. The task is to find an optimal configuration for a set of hyperbolic disks $E$ with fixed radii. We use the interior-point method as implemented in Mathematica (FindMinimum, [40]) and Matlab (fmincon, [24]).

In the most general case the problem is defined as in (3.16), where the only constraint is a geometric one, that is, the disks are not allowed to overlap. Here, the radii are fixed and the optimization concerns only the locations of the disks.

[TABLE]

This nonlinear optimization problem can be solved using the interior-point method. This solution would be a local minimum. The standard textbook reference is Nocedal and Wright [30].

Notice, that the objective function is indeed the capacity of the constellation. Often optimization problems with geometric constraints are related to packing and fitting problems. The task here is orders of magnitude more demanding since, at every point evaluation one solution of the capacity problem has to be computed, and as the disks move the constraints change as well. The number of evaluations is greater than the number of iteration steps, since the gradients and Hessians must be approximated numerically. It should be noted that the success of the optimization depends on the high accuracy of the capacity solver, since otherwise the approximate derivatives are not sufficiently accurate.

In the context of this work, there have been no attempts to devise a special method that would incorporate some of the insights gathered during this study. Instead, the numerical optimization is used to challenge those insights and therefore the optimizations have been computed with minimal input information.

4. Minimizing Capacity: Constrained Configurations

As mentioned above, even with a small number of disks the combinatorial explosion of the number of configurations is evident. Therefore, we restrict ourselves to a series of experiments each with increasing complexity building toward an understanding of the fundamental geometric principles behind the minimal configurations. In each case we consider a set of hyperbolic disks $E_{\_}j$ with radii $r_{\_}j$, where some geometric constraint is placed on all or some of the disks in the constellation.

An initial observation is that due to conformal invariance of the capacity, its numerical value remains invariant under a Möbius transformation of the unit disk onto itself. Therefore we may assume that the disk with the largest radius $r_{\_}1$ is centered at the origin.

Further, consider a disk $B_{\_}{\rho}(z_{\_}2,r_{\_}2)$ with center $z_{\_}2$ on the segment $(0,1)$. The disk lies in the lens-shaped region

[TABLE]

with $\rho_{\_}{{\mathbb{B}}^{2}}(0,\mathrm{i}v)=r_{\_}2$ where $v=\sqrt{1+\tau^{2}}-\tau$ and is tangent to both boundary arcs of $W$ and $\pm 1\in\partial W,$ see Figure 2 (right). Every disk lies within its own associated lens-shaped domain.

4.1. Disks with collinear centers

Consider a set of $m$ hyperbolic disks $E_{\_}j$ with radii $r_{\_}j$ and centers on the diameter $(-1,1)$ with $\sum_{\_}{j=1}^{m}2r_{\_}j=d_{\_}1=\rho_{\_}{\mathbb{B}^{2}}(-0.6,0.6)$. We choose the hyperbolic centers of these disks so that the hyperbolic distance between them is $d\geq 0$ where $d=0$ corresponds to the case when they touch each other. The goal is to establish upper and lower bounds for ${\rm cap}(\mathbb{B},\cup E_{\_}j)$. Since the hyperbolic radius of a hyperbolic disk is invariant under Möbius transform, in view of (2.2), we have ${\rm cap}(\mathbb{B},E_{\_}j)=2\pi/\log(1/\operatorname{\mathrm{th}}(r_{\_}j/2))$ for all $E_{\_}j$.

The cases ${\rm cap}(\mathbb{B},\cup_{\_}{j=1}^{m}E_{\_}j)$ for $m=2,3,4$ over the range $0.02\leq d\leq 4$ are shown in Figure 4. The conjectured lower bound with $d=0$ is computed with $hp$-FEM (see the ‘red dot’ in Figure 4 (right)), all other capacities are computed with BIE. From Figure 4 we also see that

[TABLE]

as the separation $d$ becomes large.

4.1.1. Verification of results

Let us consider the case with four disks and set $E=\cup_{\_}{j=1}^{4}E_{\_}j$. The initial position is when the disks are contiguous, tangent to each other, and then the hyperbolic distance $d$ between the disks increases from [math] to $0.3$. The conclusion is that the value $d=0$ yields the minimal value of the capacity of the constellation.

The values of the capacity ${\rm cap}({\mathbb{B}}^{2},E)$ in Table 2 have been computed using both methods, the FEM and the BIE method. For the BIE, we use $n=2^{7}$ and $\alpha=0.8\mathrm{i}$. Table 2 shows the absolute differences between the computed values which indicates a good agreement between the two methods. As in [13], the values computed using the FEM will be considered as reference values and used to estimate the error in the values computed by the BIE method for several values of $n$. The BIE method cannot be used for $d=0$. The error for $d=0.05,0.1,\ldots,0.3$ is presented in Figure 5 (right) which illustrates the exponential convergence with order of convergence $O(e^{-\sigma n})$ where $\sigma=-\log|\alpha|\approx 0.223$. Numerical experiments (not presented here) with other values of $\alpha$ indicate that the order of convergence depends on $\alpha$ as well as the centers $z_{\_}1,\ldots,z_{\_}m$ and the radii $r_{\_}1,\ldots,r_{\_}m$ of the inner circles. A detailed analysis of the order of convergence for the above BIE method is a subject of future work.

4.2. Four disks: Permutation of contiguous disks

We consider next two cases where all the disks of the constellation have fixed hyperbolic radii $A>B>C>D>0$ but their relative ordering is not constrained other than that each disk is tangent to at least one other disk of the constellation and their hyperbolic centers lie (a) either on the diameter $(-1,1)$ or (b) on the circle $\{z:|z|=1/2\}.$

Now the question is what is the effect of the permutation of the disks on the capacity. There are 24 permutations with 12 different capacities due to symmetry. For every realisation, the radii are denoted by $r_{\_}j$ from left to right and the constellations are denoted by $E_{\_}D$ and $E_{\_}C$, respectively. For $E_{\_}D$ we set $(A,B,C,D)=(1/2,2/5,1/4,1/5)$, and for $E_{\_}C$ slightly perturbed $(A,B,C,D)=(1/2,1/3,1/4,1/5)$. The results are collected in Table 3 and Figure 1 shows the observed extremal permutations. Interestingly, the resulting capacities have exactly the same dependence on the relative sizes of the radii.

4.3. Three immobile disks, one rolling disk

In the final experiment of the section we study the situation when one disk is free to roll on the remaining three contiguous immobile disks, centers on the diameter $(-1,1)$ and tangent to each other. The route of the mobile disk is parametrized with a parameter $\tau\in[0,1]$ where the values [math] and $1$ are for the case when also the mobile disk has its center on the diameter $(-1,1)$ and the values $1/3$ and $2/3$ correspond to the intermediate points on the route when the rolling disk is tangent to two immobile disks. Depending on the radii, it might also happen that there is only one such point. In Figure 6 below we see that for the values $1/3$ and $2/3$ the capacity of the constellation attains a local minimum. The numerical results for this example are computed using the BIE method. So, instead of assuming that the disks are touching each other, we assume that the disks are close to each other such that the hyperbolic distance between them is $d=0.02$. In all cases the hyperbolic centers of the three fixed disks are $z_{\_}1=-\operatorname{\mathrm{th}}((r_{\_}1-d)/2)$, $z_{\_}2=\operatorname{\mathrm{th}}((r_{\_}2+2d)/2)$, and $z_{\_}3=\operatorname{\mathrm{th}}((r_{\_}3+2r_{\_}2+3d)/2)$. The hyperbolic center $z_{\_}4$ of the moving disk is on the red curve shown in the figure. The observed results are summarized in the second row of Figure 6.

5. Minimizing Capacity: Optimization under Free Mobility

In this section we consider a series of experiments, where some disks are given fixed positions but the others are free to move within constraints. The constraints can restrict the admissible configurations to specific regions. In the most general case, the only constraint is that the disks should not overlap. In all simulations it is assumed that the disks have a minimal separation $\delta>0$. In those cases where the disks touch, that is, $\delta=0$, only $hp$-FEM results are reported.

5.1. Three fixed disks. One freely moving disk

Consider three hyperbolic disks with equal hyperbolic radii $=0.2$, and whose centers are at $0.5e^{2(k-1)\pi\mathrm{i}/3},\quad k=1,2,3.$ We consider a fourth hyperbolic disk whose hyperbolic radius is $r$ and its hyperbolic center is $z=x+\mathrm{i}y$ such that the four disks are non-overlapping. Let a function $u(x,y)$ be defined by $u(x,y)={\rm cap}(\mathbb{B},E),$ where $E$ is the union of the four disks. The level curves of the function $u(x,y)$ for six cases of $r$ are given in Figure 7. Notice, that the locations of the local minima depend on the chosen radius $r$ of the free disk. Due to symmetry, there is a local minimum at the origin in every case. The results suggest that there exists a critical radius $r_{\_}c$ such that the global minimum is found at the origin for all sufficiently large $r$, that is, $r>r_{\_}c$, but next to one of the fixed disks for $r<r_{\_}c$. The interior-point method is guaranteed to converge to one of the local minima, and therefore for all $r$ a local minimum may be attained when the mobile disk is centered at the origin.

5.2. One fixed disk. Two moving disks on a circle

Let us next consider three disks $D_{\_}1,D_{\_}2,D_{\_}3$ with equal hyperbolic radii $r=0.3$. The centers of these three disks are placed on the circle $|z|=0.5$. We assume that the disk $D_{\_}1$ is fixed with center on the positive real line, $D_{\_}2$ is in the upper-half plane and $D_{\_}3$ is in the lower-half plane. Starting when the three disks are touching each others (see Figure 8 (left)), these disks start moving away from each other such that the hyperbolic distance $d$ between the hyperbolic centers of $D_{\_}1$ and $D_{\_}2$ is the same as for $D_{\_}1$ and $D_{\_}3$. When all these disks are touching each other, $d=2r$. The maximum value $d_{\_}{\rm max}$ of $d$ is obtained when the the disks $D_{\_}2$ and $D_{\_}3$ are touching each other (see Figure 8 (middle)). The values of the capacity as a function of $d$ are shown in Figure 8 (right) where the values of the capacity for $2r<d<d_{\_}{\rm max}$ are computed by the BIE method and; for $d=2r$ and $d=d_{\_}{\rm max}$ by the FEM. The minimal capacity is found when $d=2r$ and the maximal when the centers of the three disks form an equilateral triangle.

5.3. One fixed disk. Three moving disks on a circle

Staying on the circle $|z|=0.5$ we consider four disks with centers on the circle and hyperbolic radii $3/30$, $5/30$, $7/30$, and $9/30$. Without any loss of generality, we will assume that the disk with hyperbolic radius $9/30$ is fixed with its center on the positive real line at the point $0.5$. Then, we search for the positions of the other three disks that minimize the capacity. The initial positions of these three disks are assumed to be $0.5e^{2k\pi\mathrm{i}/4}$ for $k=1,2,3$. For the optimized positions, we have obtained six positions, with three different values of the capacity due to symmetry (see Figure 9). For the disks in the first column in Figure 9, the capacity is $4.6269$. The capacity is $4.6193$ for the second column and $4.6621$ for the third column.

5.4. One fixed disk. Three moving disks

Finally, we consider four disks with hyperbolic radii $3/30$, $5/30$, $7/30$, and $9/30$. This time, we will assume that the disk with hyperbolic radius $9/30$ is fixed with its center at the origin. The task is to find the positions of the three free disks that minimize the capacity where the initial positions of the three disks are assumed to be $0.5e^{2k\pi\mathrm{i}/3}$ for $k=0,1,2$. For the optimized positions, we have obtained two configurations, as shown in Figure 10, with the capacity $4.2322$ which is the global minimum.

If we assume that the three disks with hyperbolic radii $5/30$, $7/30$, and $9/30$ have fixed positions as in Figure 10 (left), and the small disk with hyperbolic radius $3/30$ is moving. Assume that the center of the small disk is $z=x+\mathrm{i}y$ such that the four disks are non-overlapping. Let a function $u(x,y)$ be defined by $u(x,y)={\rm cap}(\mathbb{B},E),$ where $E$ is the union of the four disks. The level curves of the function $u(x,y)$ are given in Figure 10 (right). As we can see from the figure, the capacity has three local minima and the capacity for the position in Figure 10 (left) is the global minimum. This experiment has been repeated multiple times with different initial starting positions for the free disks and every one one of the local minima has been observed.

5.4.1. On Computational Costs

Naturally, the optimisation problems are the computationally most expensive ones of all our numerical experiments. In Table 5 performance data on the four disks free mobility problem is presented. Comparison of the two methods is only qualitative, since both underlying hardware and the interior-point implementations are different. However, some conclusions can be derived. In all cases the interior-point tolerance is the same, $\epsilon=10^{-6}$, and within the $hp$-FEM simulations, meshing is performed with the same discretization control in every evaluation. For optimal performance, the individual solutions must be accurate enough so that the error induced by numerical approximation of the gradients and Hessians is balanced with other sources of error. For the $hp$-FEM it appears that the same mesh with $p=4$ is not adequate in comparison with the one at $p=8$. Even though the time spent in one individual iteration step is doubled, the overall time for $p=8$ is significantly lower. At every evaluation the number of degrees of freedom is roughly 13000 (initial configuration: 13542, and final: 12589). Similarly, for BIE the performance at $n=2^{7}$ is superior to that at $n=2^{4}$.

The two implementations have very different requirements per iteration step. Observe that the number of iteration steps is comparable, yet the number of evaluations is not. The average time for one evaluation in BIE is four to five times faster than one evaluation in $hp$-FEM. Matlab and Mathematica results have been computed on modern Intel and Apple Silicon computers, respectively.

5.5. Hyperbolic area lower bound

Finally, we compute the capacity of a constellation of disjoint hyperbolic disks and compare the computed values with the Hyperbolic area lower bound [9]. Let $E_{\_}r$ be the union of $m$ disjoint hyperbolic disks with equal hyperbolic radii $r$ such the hyperbolic distance between any two disks is $0.02$ (see Figure 11 for $r=0.5$ and $m=2,3,4$). For $m=4$, we consider two cases (as shown in Figure 11) where the centers of the disks in Case I are on the real and imaginary axes. In Case II, the centers are on the rays $e^{\mathrm{i}\theta}$ for $\theta=0,\pi/3,2\pi/3,4\pi/3$. The hyperbolic area of these $m$ disks is $4m\pi\operatorname{\mathrm{sh}}^{2}(r/2)$. Consider the hyperbolic disk $B_{\_}\rho(0,M)$ whose hyperbolic area is the same as the hyperbolic area of the $m$ disks, then

[TABLE]

Then $L(r)=\operatorname{\mathrm{cap}}(\mathbb{B},B_{\_}\rho(0,M))$ is the hyperbolic area lower bound of $\operatorname{\mathrm{cap}}(\mathbb{B},E_{\_}r)$. In view of (2.3), we have

[TABLE]

The BIE method is then used to compute $\operatorname{\mathrm{cap}}(\mathbb{B},E_{\_}r)$ for several values of $r$ with $0.02\leq r\leq 2$. Our computed minimum value of the capacity can be considered a lower bound of the capacity of the constellation of $m$ disjoint hyperbolic disks. We compare the computed value with the hyperbolic area lower bound by defining

[TABLE]

The graph of $L_{\_}r$ is shown in Figure 12 for $0.02\leq r\leq 2$ and $m=2,3,4$. As $r\to\infty$ it appears that the improvement tends to zero. This is a consequence of the nature of hyperbolic geometry. With one disk fixed in the centre the other three will have ever smaller contributions to the capacity since their Euclidean areas tend to zero as in Figure 2 (right). It is an indication of the complexity of the problem that the graphs in Figure 12 do not reveal any simple connection between the number of the disks and the minimal capacity.

6. Conclusions

We study lower bounds for the conformal capacity of a constellation of disjoint hyperbolic disks $E_{\_}j\subset\mathbb{B}^{2}$, $j=1,...,m$, using a novel idea: instead of using a symmetrization transformation, which usually leads to fusion of the disjoint disks, we are looking for a lower bound in terms of another constellation which yields a minimal value. The traditional symmetrization transformation [31], [3], [8], is now replaced by free mobility of individual disks with the constraint that the hyperbolic radii of the disks are invariant and the disks are non-overlapping. In this process, due to the conformal invariance, the conformal capacity of each disk stays invariant, whereas the capacity of the whole constellation may significantly vary. Moreover, the hyperbolic area of the constellation is also constant.

The optimization methods we used produced (locally) minimal constellations such that the disks group together, as closely as possible. This coalescing is reminiscent of the behavior of some animal colonies in cold weather conditions for the purpose of heat flow minimization. Mathematical methods are not available for analytic treatment of the problems, but we are convinced that there is a strong connection with combinatorial geometry, topics like packing and covering problems. Such problems often have many local minima [6, p. 157].

We carried out numerical simulations using two different methods, the BIE and $hp$-FEM methods and the close agreement of the two computational methods confirmed the results. Because of the complexity of the problem we studied various subproblems where disk centers satisfied constraints such that the centers are on the interval $(-1,1)$ or at the same distance from the origin. In both cases we observed the grouping phenomenon (cf. Figure 1) and, moreover, noticed that permutation of disks has influence on the capacity if the radii are different. Because the hyperbolic area of a constellation is a constant, it is now clear that the hyperbolic area alone does not define the constellation capacity.

This observation led us to compare our computed lower bound to Gehring’s sharp lower bound given in terms of hyperbolic area. The conclusion was that we obtained in some cases approximately $10\%$ improvement when $m=4$.

The numerical agreement of the BIE and $hp$-FEM methods was very good, typically ten decimal places or better, and the expected exponential convergence was observed, see Figure 5. The performance of the BIE method was significantly faster than the $hp$-FEM method when it comes to computational time and flexilibity to modify the code to new situations. This is probably due to the heavy data structure of the $hp$-FEM method due to hierarchial triangulation refinement process of the method.

A vast territory of open problems remains. First, it would be interesting to study whether some kind heuristic methods would lead to "close to extremal" constellations, to be used as initial steps of the minimization. Such a method could be based on some computationally cheaper object function than the capacity itself: for instance, first, the maximization of the number of the mutual contact points of the constellation. Second, the case of $m>5$ disks of equal radii seems to be completely open. Perhaps in this case the number of locally minimal constellations grows exponentially as a function of $m.$ Third, one could study constellations of other types of geometric figures like hyperbolic triangles.