Computation of conformal invariants

We study numerical computation of conformal invariants of domains in the complex plane. In particular, we provide an algorithm for computing the conformal capacity of a condenser. The algorithm applies for wide kind of geometries: domains are assumed to have smooth or piecewise smooth boundaries. The method we use is based on the boundary integral equation method developed and implemented in [N3]. A characteristic feature of this method is that, with small changes in the code, a wide spectrum of problems can be treated and we include code snippets within the text to indicate implementation details. We compare the performance and accuracy to previous results in the cases when numerical data is available and also in the case of several model problems where exact results are available.


Introduction
During the past fifty years conformal invariants have become crucial tools for complex analysis. Most important of these invariants are the harmonic measure, the conformal capacity, the extremal length, and the hyperbolic distance [Ah, GM, KL, ST]. But this is not all: the generalized capacity, the transfinite diameter, the reduced extremal length, the hyperbolic area, and the modulus metric [Du,Ki,Va,V1,V2] are some additional examples. Some of the many applications of these tools are discussed in the articles of the handbook [Ku]. In view of the plenitude of these applications, it is surprising that these invariants can be expressed explicitly only in very few special cases. Sometimes rudimentary upper or lower bounds for conformal invariants in terms of less involved comparison functions can provide important steps in proofs.
At the same time it seems that the full power of conformal invariance remains unused. One reason for this is that the analytic expressions for conformal invariants are usually too complicated for pen and paper calculations and the existing computational methods are scattered throughout the mathematical literature: the way from theory to practical experimentation is too long. On the other hand, the creators of the existing computational methods may not be aware of the scope of applicability of their methods in theoretical studies: the way from experiments to theory is also long. If the distance from theory to experiments could be made shorter, a theoretical researcher could easily experiment the dependence of a problem under perturbation of geometry and vice versa a computational scientist would find new types of benchmark problems and areas of application.
The above ideas as our guiding principles, we have written a series of papers of which this is the first one [KNV, NV]. As far as we know, our work is the first attempt to provide computational tools for a wide class of conformal invariants with the feature that modification of geometry is simple. The method we use was developed and implemented by the first author [N3] and we apply it to study several computational problems never studied before and we also compare its performance to several results in the literature. As test problems we use the computation of condenser capacity, a topic studied by the second author in several papers [HRV1,HRV2,HRV3].
A condenser is a pair (D, C) where D is an open set in R 2 and C ⊂ D is compact. In our study we assume that the topology is simple but still general enough for most applications: the sets ∂D and ∂C are connected sets, each set is a piecewise smooth Jordan curve.
The conformal capacity, or capacity for short, is defined by [Du] (1.1) cap(D, C) = D |∇u| 2 dm where u : D → R is a harmonic function with u(x) ≥ 1 for all x ∈ C and u(x) → 0 for x → ∂D . The domain G = D \ C is called the field of the condenser and the closed sets C and D c are called the plates of the condenser. Then, the capacity cap(D, C) may alternatively be written as cap(G).
In literature, only very few formulas are given for the capacities of concrete condensers. Numerical methods are therefore needed to compute the value of (1.1). Our problem is reduced to the classical problem of solving numerically the Dirichlet boundary value problem for the Laplace equation. Moreover, by the Dirichlet principle [GM,, the extremal function u 0 is harmonic and minimizes the integral [Ah], [GM,: where the infimum is taken over all C o (D) functions with the indicated Dirichlet boundary values. The capacity of condensers is invariant under conformal mappings, and hence domains with difficult geometry can be handled with the help of conformal mappings [AVV,DT,Du,PS,SL,Va,V2]. See also [BBG, DEK, ET].
Before proceeding to the contents of our work a few general remarks about the literature we know about may be in order. Because of the wide scope of conformal invariants, relatively few cases exist where "the right answer" is known. In such cases, the computational performance may be analysed by observing convergence features of the results under successive refinements of the numerical model, and error estimates maybe based on general theory. In those relatively few cases we have found in the literature where the analytic formula is known, the true error estimate may be given. Sometimes a high accuracy can be achieved, say 12 decimal places, but the dilemma is that if the geometry of the problem is smoothly changed a bit, the method might not be applicable at all.
Section 2 summarizes our computational workhorse, the boundary integral method geared for the capacity computation of ring domains, which will be applied in several later sections, sometimes together with auxiliary procedures. In Sections 3 and 5 we consider ring domains for which the exact value of capacity is known and investigate the performance of our method. Sections 4, 6, and 7 deal with condensers whose one or two complementary components are slits-these are well-known examples computationally challenging problems and we use here auxiliary conformal mapping to overcome computational difficulties. Section 8 deals with the case when both complementary components of a ring domain are thin rectangles. In Section 9, we consider the numerical computation of the hyperbolic capacity and the elliptic capacity of compact and closed sets. The final Section 10 gives some concluding remarks and information about the access to our MATLAB software. A domain G in the extended complex plane C = C ∪ {∞} , whose complement C \ G has two components, is called a ring domain or, briefly, a ring. It is a classical fact that a ring can be mapped by a conformal map onto an annulus {z : q < |z| < 1} , q ∈ (0, 1) . A ring R is the simplest example of a condenser and its capacity is given by [Du], [GM, cap(R) = 2π log(1/q) .
Because of the conformal invariance of the capacity, this definition is independent of the conformal map. For the computation of the capacity we will often use an auxiliary conformal mapping to avoid computational singularities.
In this section we describe the method of our numerical work, based on the solution of the boundary integral equation with the generalized Neumann kernel [N3,WN]. The integral equation has been applied to calculate conformal mappings onto several canonical domains [N1,N4,NF]. We review the application of the integral equation to compute the conformal mapping from doubly connected domains onto an annulus {z : q < |z| < 1}, q ∈ (0, 1), and present the MATLAB implementation of the method. In later sections we will apply this method for capacity computation of several condensers, in particular we will consider several types of rings with a simple geometric structure.
2.3. The generalized Neumann kernel. Let G be a bounded or an unbounded doubly connected domain bordered by Γ = ∂G = Γ 1 ∪ Γ 2 where each of the boundary components Γ 1 and Γ 2 is a closed smooth Jordan curve. We choose the orientation of boundary Γ such that when we proceed along Γ , the domain G is always on its left. If G is bounded, then Γ 1 is the external boundary and Γ 2 is contained in the bounded domain whose boundary is Γ 1 . The complement C\G of the domain G with respect to the extended complex plane C consists of two simply connected domains G 1 on the right of Γ 1 and G 2 on the right of Γ 2 . The domain G 2 is bounded and the domain G 1 is unbounded with ∞ ∈ G 1 . Further, we assume α is an auxiliary given point in the domain G and z 2 is an auxiliary given point in the simply connected domain G 2 . When G is unbounded, then ∞ ∈ G and the two domains G 1 and G 2 are bounded. We assume that z 1 and z 2 are auxiliary given points in the simply connected domains G 1 and G 2 , respectively. See Figure 1. We parametrize each boundary component Γ j by a 2π-periodic complex function η j (t), t ∈ J j := [0, 2π], j = 1, 2. We assume that each of these functions η j (t) is twice continuously differentiable with η ′ (t) = 0 (the presented method can be applied also if the curve Γ j have a finite number of corner points but no cusps [NMZ]). Then we define the total parameter domain J as the disjoint union of the two intervals J 1 = [0, 2π] and J 2 = [0, 2π], i.e., The elements of the total parameter domain J are ordered pairs (t, j) where t is a real number, 0 ≤ t ≤ 2π, and the index j is an integer indicates that the interval J j containing t [N3]. Hence, the boundary Γ can be parametrized by (2.4) η(t, j) = η j (t), t ∈ J j , j = 1, 2.
For a given t, the index k such that t ∈ J k will be always clear from the context, see e.g., [N1,N3,NF,NG,WN]. So the pair (t, k) in the left-hand side of (2.4) will be replaced by t and a parametrization of the whole boundary Γ can be defined on J by We denote by H to the space of all functions of the form where ρ 1 (t) and ρ 2 (t) are 2π-periodic Hölder continuous real functions on J 1 and J 2 , respectively. Let A be the complex function [N3] (2.6) where θ is a real function with constant value on each interval J j , i.e., The kernel N(s, t) is continuous [WN]. Hence, the integral operator N defined on by is compact. The integral equation with the generalized Neumann kernel involves also the following kernel which has a singularity of cotangent type [WN]. The integral operator M defined on H by is singular, but is bounded on H [WN]. For more details, see [WN]. For the above function A defined by (2.6), the following integral equation is uniquely solvable for any real function γ in H [N2]. Furthermore, if ρ is the unique solution of the boundary integral equation (2.9), then the real function h defined by is a piecewise constant function on the boundary Γ, i.e., where h j is a real constant, j = 1, 2 [N2]. Moreover, are boundary values of an analytic function f in the domain G with f (∞) = 0 for unbounded G. For more details, see [N2, N3] and the references cited therein.
2.12. Numerical solution of the boundary integral equation. The MATLAB function fbie in [N3] provides us with a fast and efficient method for solving the boundary integral equation (2.9). The function fbie is based on discretizing the boundary integral equation (2.9) using the Nyström method with the trapezoidal rule [AKT, At, TW]. This discretization leads to a non-symmetric linear system. Then, the MATLAB function gmres is used to solve the linear system. The matrix-vector multiplication in the GMRES method is computed using the MATLAB function zfmm2dpart in the toolbox FMMLIB2D [GG]. The function fbie provides us also with approximations to the piecewise constant function h in (2.10). The computational cost for the overall method is O(n log n) operations where n (an even positive integer) is the number of nodes in each of the intervals J 1 and J 2 . For the accuracy of the obtained numerical results, it is known that the order of the convergence of the Nyström method depends on the order of convergence of the used quadratic method [At]. The quadratic method used in the function fbie is the trapezoidal rule which gives surprisingly accurate numerical results for periodic functions [At, TW]. In view of of (2.7) and (2.8), the smoothness of the two kernels N(s, t) and M(s, t) depends on the smoothness of the parametrization function η(t). Similarly, in this paper, the function γ on the right-hand side of the integral equation (2.9) will be defined in terms of the parametrization η(t). Hence, the smoothness of the function γ will depends also on the smoothness of the boundary Γ. Thus, the order of convergence the trapezoidal rule depends on the smoothness of the boundary Γ of the domain G. For domain with smooth boundaries, we use the trapezoidal rule with equidistant nodes. Thus, for domains with C ∞ smooth boundaries, the integrand in the integral equation (2.9) will be C ∞ smooth and hence the rate of convergence of the numerical method is O(e −cn ) with a positive constant c (see [K2, p. 223]). If the boundary is C q+2 smooth (q ≥ 0), then the function γ will be also C q+2 smooth and the rate of convergence of the numerical method is O(1/n q ) [K1]. 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 [K1]. The graded mesh is obtained by substituting a new variable to remove the discontinuity in the derivatives of the solution ρ(t) of the boundary integral equation (2.9) at the corner points [K1,LSN].
To use the MATLAB function fbie, the vectors et, etp, A, and gam that contain the discretizations of the functions η(t), η ′ (t), A(t), and γ(t), respectively, will be stored in MATLAB. Then we call the function [rho, h] = fbie (et, etp, A, gam, n, iprec, restart, gmrestol, maxit) to compute the vectors rho and h which contain the discretizations of the solution of the integral equation ρ(t) and the piecewise constant function h(t), respectively. In the numerical experiments in this paper, we set the tolerances of the FMM and the GMRES to be 0.5 × 10 −15 and 10 −14 by choosing iprec = 5 and gmrestol = 10 −14 , respectively. We use the GMRES without restart by choosing restart = [ ]. The maximum number of iterations for GMRES is chosen to be maxit = 100.
Choosing the value of n depends on the geometry of the domain G. If the domain G has a simple geometry and smooth boundary, we can obtain accurate numerical results by choosing moderate values of n. If the domain G has a complex geometry such that its boundary has corners or its boundary components are close to each other, it is required to choose a sufficiently large value of n to obtain accurate results. For domains with corners, we choose n as a multiple of the number of corners.
Once the discretizations of the two functions ρ(t) and h(t) are computed, we can compute a discretization of the boundary values of the analytic function f (z) through Then approximations to the values of the function f (z) for any vector of points z in G can be obtained using the Cauchy integral formula. Numerically we carry out this computation by applying the MATLAB function fcau [N3] by calling fz = fcau (et, etp, fet, z) for bounded G and by calling fz = fcau (et, etp, fet, z, n, 0) for unbounded G (here 0 = f (∞)).
For more details, we refer the reader to [N3]. The computations in this paper 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 MATLAB tic toc commands were used to measure the computation times.
2.13. Computing the conformal mapping for bounded domains. If the domain G is bounded, then can compute the conformal mapping w = Φ(z) from G onto the annulus {w ∈ C : q < |w| < 1} with the normalization Φ(α) > 0 as in the following theorem from [N1]. Here, α is an auxiliary given point in G.
Theorem 2.14. Let θ 1 = θ 2 = π/2, let the function A be defined by (2.6), and let the function γ be defined by If ρ is the unique solution of the boundary integral equation (2.9) and the piecewise constant function h is given by (2.10), then the function f with the boundary values (2.11) is analytic in the domain G, the conformal mapping Φ is given by and the modulus q is given by 2.18. Computing the conformal mapping for unbounded domains. For unbounded domain G, the following theorem from [N1] provides us with a method to compute the conformal mapping w = Φ(z) from G onto the annulus {w ∈ C : q < |w| < 1} with the normalization Φ(∞) > 0.
Theorem 2.19. Let θ 1 = θ 2 = π/2, let the function A be defined by (2.6), and let the function γ be defined by If ρ is the unique solution of the boundary integral equation (2.9) and the piecewise constant function h is given by (2.10), then the function f with the boundary values (2.11) is analytic in the domain G with f (∞) = 0, the conformal mapping Φ is given by and the modulus q is given by 2.23. Computing the capacity of the doubly connected domain G . Since the capacity is invariant under conformal mapping, we shall compute the capacity of the above doubly connected domain G (for both cases, bounded and unbounded) by mapping G onto the annulus R = {w ∈ C : q < |w| < 1} using the method presented in the above two theorems. Then the capacity of G is the same as the capacity of the annulus R which is given by The method presented in this section for computing the radius q of the inner circle of the annulus R = {w ∈ C : q < |w| < 1} and hence the capacity cap(G) for both bounded and unbounded doubly connected domains G can be implemented in MATLAB as in the following function. function [q,cap] = annq (et,etp,n,zz,z2,type) % This function computes the inner radius q for the conformal mapping % w=Phi(z) from doubly connected domain G onto the annulus q<|w|<1 and % the capacity of G, cap(G)=2pi/log(1/q), where: % et, etp: the parametrization of the boundary of G and its derivative % n: the number of discretization points % zz: zz=alpha is a given point in $G$ for bounded G; % and zz=z1 is a given point in the interior of the curve that will be % maped onto the unit circle for unbounded G % z2: a given point interior to the curve that will be maped onto the % circle |w|=q for both cases of bounded and unbounded G. % type='b' for bounded G; and type='u' for unbounded G % if type=='b' alpha = zz; A = et-alpha; gam = -log(abs ((et-z2)./(alpha-z2))); elseif type=='u' z1 = zz; A = ones(size(et)); gam = -log(abs ((et-z2)./(et-z1))); end [~,h] = fbie (et,etp,A,gam,n,5,[],1e-14,200); q = exp(mean(h(n+1:2 * n))-mean(h(1:n))); cap = 2 * pi/log(1/q); end

Rings with piecewise smooth boundaries
In this section, we apply the method presented in Section 2 to compute the capacity of several doubly connected domains G with piecewise smooth boundaries. For the first two examples, the exact value of the capacity is known.
3.1. Two circles. In th first example, we consider the doubly connected domain in the exterior of the two circles |z| = 1 and |z − a| = r where r > 0 and a is a real number with a > 1 + r.
Using a Möbius transform, the domain G can be mapped onto the annulus {w ∈ C : q < |w| < 1} where q is obtained by solving the following equation which follows from the invariance of the cross ratios under Möbius transformations [V2, (1.28)] Hence, the exact value of conformal capacity is given by cap(G) = 2π/ log(1/q) . We use the MATLAB function annq with z 1 = 0, z 2 = a and n = 2 10 to compute approximate values for the capacity for several values of a and r. First, we fixed r = 1 and chose values of a between 2.01 and 6. The relative errors for the computed values 3.2. Two confocal ellipses. In this example, we consider the bounded doubly connected domain G in the interior of the ellipse and in the exterior of the ellipse where r 1 > r 2 > 1. The domain G is the image of the ring q = r 2 /r 1 < |ζ| < 1 under the Joukowski map Hence, the exact value of conformal capacity of G is given by cap(G) = 2π/ log(1/q) = 2π/ log(r 1 /r 2 ). We use the MATLAB function annq with α = ((r 1 + 1/r 1 ) + (r 2 + 1/r 2 ))/4 ∈ G, z 2 = 0, and n = 2 12 to compute approximate values for the capacity for several values of r 1 and r 2 . First, we fixed r 2 = 2 and chose values of r 1 between 2.05 and 6. The relative errors for the computed values are presented in Figure 3 (left). Then, we fixed r 1 = 4 and chose values of r 2 between 1.01 and 3.9. The relative errors for the computed values for this case are presented in Figure 3 (right). 3.3. Complete elliptic integrals. We recall the following facts about complete elliptic integrals and hypergeometric functions, needed for the sequel. The Gaussian hypergeometric function is the analytic continuation to the slit plane C \ [1, ∞) of the series where a, b, and c are complex numbers with c = 0, −1, −2, . . .. Here (a, 0) = 1 for a = 0, and (a, n) is the Appell symbol or the shifted factorial function (a, n) = a(a + 1)(a + 2) · · · (a + n − 1) The complete elliptic integrals, K(r) and K ′ (r), of the first kind are The elliptic integrals, E(r) and E ′ (r), of the second kind are Then K : (0, 1) → (0, ∞) is an increasing homeomorphism and E : (0, 1) → (1, π/2) is a decreasing homeomorphism. The decreasing homeomorphism µ : The basic properties of these functions can be found in [AVV, BF, OLBC]. For example, it follows from [AVV,(5.2)] for 0 < r < 1 that In the numerical calculations in this paper, we compute the values of µ(r) through (3.7) where the values of K(r) and K ′ (r) are computed by the MATLAB function ellipke.
Since 0 < r < 1 and r ′ = √ 1 − r 2 , we can show that r < 2 √ r 1 + r < 1 and 0 < 1 − r ′ 1 + r ′ < r. Thus, when r is too close to 0, we can use the first formula in (3.8) to get accurate results with MATLAB function ellipke. When r is very close to 1, we use the second formula in (3.8).
3.9. Jacobi's inversion formula for µ. In his fundamental work on elliptic functions, C.G.J. Jacobi proved several dozens of formulas for these functions and related functions involving infinite products or theta functions. As pointed out in [AVV,Thm 5.24(2)], some of these formulas can be rewritten so as to give formulas for µ −1 (y) . We give two examples. Jacobi's inversion formula for µ is [AVV,Thm 5.24 (2)] Another example of Jacobi's work is the following formula for µ −1 (y) in terms of theta functions Because these theta functions converge very fast in [0, 0.95] , a few terms of series expansion are enough to achieve numerical values correct up to 15 decimal places. A Newton algorithm for computing µ −1 (y) was implemented in [AVV,pp. 92,438].
3.12. Square in square. In our third example, we consider the doubly connected domain G which is the difference of two concentric squares where 0 < a < 1.
The exact value of the capacity of this domain is given by [Bo, (3.13) cap(G) = 4π µ(r) , Then, by [AVV,Exercises 5.8(3)], we have By [AVV,(5.2)], we have Thus, it follows from (3.13) that (3.14) cap(G) = 8 π µ(2uv) , We use the MATLAB function annq with α = 1 + a ∈ G, z 2 = 0, and n = 2 17 to compute approximate values for the capacity for several values of a between 0.1 and 0.9. The obtained results are presented in Table 1. Table 1 presents also the exact values of the capacity and the numerical results computed in [HRV1] using an hp-FEM algorithm. We see from the results presented in the table that accurate results can be obtained using the presented method. The last column in Table 1 presents the CPU time (in seconds) for our method. The GMRES requires between 23 to 25 iterations only to converge. The obtained results using the presented method are not as accurate as the results obtained by the hp-FEM algorithm in [HRV1]. This is expected when we compare BIM and FEM for domains with corners.
3.15. Polygon in polygon. In the forth example, we consider the doubly connected domain G between two polygons. We assume that both polygons have m vertices where  2.5 5 9.62720096044514 9.6266 2.6 7 9.25977557690559 9.2598 2.4 9 9.15441235751744 9.1541 2.1 15 9.08360686195382 1.8 30 9.06705650051687 1.5 m ≥ 3. We assume that the vertices of the external polygon are the roots of the unity and hence lie on the unit circle |z| = 1. For the inner polygon, we assume that the vertices are the roots of the unity multiplied by q = 0.5 and thus lie on the circle |z| = q (see Figure 4 for m = 5).
The exact value of capacity of the domain is unknown (except for m = 4 where the capacity can be computed as in the square in square example, which for q = 0.5, is 10.2340925693681). We use the MATLAB function annq with α = (1 + q)/2 ∈ G, z 2 = 0, and n = 40320 to compute approximate values for the capacity for several values of m. The computed capacity is presented in Table 2. As we can see from the table, as m increases, the capacity approaches the capacity of the annulus q < |z| < 1 which is 2π/ log(1/q). For q = 0.5, the capacity of the annulus is 9.064720283654388. For some values of m, Table 2 presents also approximate values of the capacity from [BSV]. The last column in Table 2 presents the CPU time (in seconds) for our method.

Complement of two slits
In this section, we consider the doubly connected domain Ω whose complementary components are the two nonintersecting segments E = [a, b] and F = [c, d] where a, b, c and d are complex numbers (see Figure 5 (left) for a = 0, b = 1, c = 1 − i and d = 3 + 2i). Computing the capacity of such domain Ω has been considered recently in [DNV] using Weierstrass elliptic functions. Here, we shall compute the capacity of Ω using the method presented in Section 2. However, a direct application of the method presented in Section 2 is not possible since the boundaries of Ω are not Jordan curves. So, we need to first map this domain Ω onto a domain G of the forms considered in Section 2. Up to the best of our knowledge, there is no analytic formula for a conformal mapping from the above doubly connected domain Ω onto a doubly connected domain G bordered by smooth Jordan curves. So, we need to use numerical methods to find such an equivalent domain G. An iterative method for computing such a domain G has been presented recently in [NG]. Using this iterative method, a conformally equivalent domain G bordered by ellipses can be obtained as in Figure 5 (right). For details on this iterative numerical method for computing the domain G, we refer the reader to [NG]. For the new domain G, we use the MATLAB function annq with n = 2 11 to compute the capacity of G for several values of the constants a, b, c and d, as in the following examples.
4.1. Two segments on the real axis. When E = [0, 1] and F = [c, d] with d > c > 1 are real numbers, the exact capacity of Ω is known and is given by [V2, 5.54 (1), 5.60(1)] . We tested our methods for several values of c and d. First, we fixed c = 2 and chose d between 2.1 and 10. Then we fixed c = d − 1 and chose d between 2.1 and 10. The relative errors of the computed values for this case are presented in Figure 6. As we can see from Figure 6, the presented method gives accurate results with relative error around 10 −14 . Table 3 presents the approximate values of the capacity, the exact values of the capacity, and the total CPU time for several values of c and d.

Rings with a segment as a boundary component
In this section, we shall compute the capacity of doubly connected domains Ω whose boundary components are a slit and a piecewise smooth Jordan curve. Such domains cannot be mapped directly onto an annulus using the method presented in Section 2. To use the method presented in Section 2, we shall use first elementary mappings to map the domain Ω onto a domain G of the types considered in Section 2. Then the domain G is mapped onto an annulus R = {z ∈ C : q < |z| < 1} and hence the capacity of Ω is    2π/ log(1/q). In this subsection we consider two examples where the exact value of the capacity for the first example is known. 5.1. Segment and circle. In this example, we consider the doubly connected domain Ω in the exterior of the segment [0, 1] and the circle |z − a| = r where r > 0 and a is a real number with a > 1 + r (see Figure 12(left) for a = 2 and r = 0.9).  The exact value of conformal capacity of the given domain Ω is known and given by [V2, 5.54 (2)] Figure 12. The segment and circle domain Ω for a = 2 and r = 0.9 (left); and the image of this unbounded domain under the mapping ζ = Ψ −1 (z) (right). and µ is given by (3.7).
To apply our method presented in Section 2, we shall use first elementary mappings to map the domain Ω onto a domain G of the types considered in Section 2. It is known that the function z = Ψ(ζ) = 1 4 ζ + 1 ζ + 1 2 maps conformally the interior of the unit circle |ζ| = 1 onto the exterior of the segment [0, 1]. Hence, its inverse function is given by where we choose the branch for which √ 1 = 1. The function ζ = Ψ −1 (z) maps the segment [0, 1] onto the unit circle |ζ| = 1 and the exterior of the segment [0, 1] onto the interior of the unit circle |ζ| = 1. The function ζ = Ψ −1 (z) maps also the circle |z − a| = r in the z-plane onto a smooth Jordan curve inside the unit circle |ζ| = 1. Consequently, the function ζ = Ψ −1 (z) maps the doubly connected domain Ω onto a bounded doubly connected domain G of the form considered in Section 2 (see Figure 12 (right)).
Then for the domain G, we use the MATLAB function annq with n = 2 11 to compute approximate values for the capacity of Ω for several values of a and r. First, we fixed r = 1 and chose values of a between 2.05 and 6. The relative errors of the computed values are presented in Figure 13 (left). Then, we fixed a = 4 and chose values of r between 0.05 and 2.95. The relative errors of the computed values for this case are presented in Figure 13 (right). We also presents the approximate values and the exact values of the capacity for several values of r and a in Table 6.
and 1 < c < d (see Figure 14 (center)). For r = 0, F r reduced to the segment F 0 = [c, d] and hence G 0 is the doubly connected domain exterior to the two segments E = [0, 1] and F 0 = [c, d] as considered in Subsection 4.1 (see Figure 14 (left)). The exact value of cap(E, F 0 ) is given by (4.2), i.e., where µ is given by (3.7). When r = b, F b is the closed disk |z − a| ≤ b and the domain G b is then the doubly connected domain exterior to the segment E = [0, 1] and the closed disk F b = {x | |z − a| ≤ Figure 14. The domain G r with c = 1.5, d = 3.5 for r = 0 (left), r = 0.1 (center) and r = b (right).
b} which considered in Subsection 5.1 (see Figure 14 (right)). The exact value of cap(E, F b ) is given by (5.2), i.e., It is clear from the definition of the closed set F r that F 0 ⊆ F r ⊆ F b for 0 ≤ r ≤ b. As r changes continuously from 0 to b, the closed set F r changes continuously from the segment F 0 to the disk F b . Here, we shall compute the exact value of the capacity cap(E, F r ) and show that cap(E, F r ) will change from cap(E, F 0 ) to cap(E, F b ) as r changes from 0 to b.
By the elementary mapping the unbounded domain G is mapped conformally onto the unbounded domain G 1 exterior to the segment [−a/b, −(a − 1)/b] and the ellipse We can easily show that the function ζ = Ψ 2 (ξ) = ξ + 1 − (r/b) 2 4 1 ξ maps the domain exterior to the circle |ξ| = (1 + r/b)/2 onto the domain exterior of the ellipse. Hence, the inverse mapping maps the domain G 1 onto the domain G 2 exterior to the circle |ξ| = (1 + r/b)/2 and the segment [c 1 , d 1 ] where (5.7) and the branch of the square root is chosen such that √ 1 = 1. Finally, the function maps the domain G 2 onto the domain G 3 exterior to the circle |w −â| =r and the segment [0, 1] where Hence, the analytic value of cap(E, F r ) can be obtained since the exact value of conformal capacity of the domain G 3 is known [V2,5.54(2)], (5.9) cap(E, F r ) = 2π µ(τ r ) where (5.10) τ r =r a 2 −â −r 2 . The value of τ r can be obtained in terms of c, d and r as following For r = 0, the capacity given by (5.9) becomes After tedious algebra, we find that s in (5.5) is related to τ 0 in (5.12) through which, in view of (3.8), implies that µ(τ 0 ) = 2µ(s).

Hence,
cap(E, F 0 ) = π µ(s) = 2π µ(τ 0 ) and thus the capacity cap(E, F r ) given by (5.9) reduced to the capacity cap(E, F 0 ) for r = 0. Furthermore, when r = b, then it follows from (5.7) and (5.8) thatâ = a and r = r = b. Hence, it follows from (5.10) that  which implies that the capacity cap(E, F r ) given by (5.9) reduced to the Formula (5.6) for r = b. The values of cap(E, F r ) for c = 1.5, d = 3.5, 0 ≤ r ≤ b (where b = (d−c)/2 = 1) is given in Figure 15. As we can see from Figure 15, the capacity cap(E, F r ) changes continuously and rapidly increases from cap(E, F 0 ) to cap(E, F b ) as r changes continuously from 0 to b. 5.14. Segment and polygon. In this example, we consider the doubly connected domain Ω in the exterior of the segment [0, 1] and a polygon with m vertices where m ≥ 3. We assume that the vertices of the polygon are given by v k = a − re −2kπi m , k = 0, 1, 2, . . . , m − 1 (see Figure 16 (left) for a = 1.6, r = 0.5, and m = 3).
The exact value of conformal capacity is unknown. To use the method presented in Section 2, we first use the mapping function ζ = Ψ −1 (z) in Subsection 5.1 to map the doubly connected domain Ω in the exterior of the segment [0, 1] and the polygon onto a r a m Capacity (segment and polygon) Capacity (segment and circle)

The upper half-plane with a slit
In this section, we consider the doubly connected domain Ω = H 2 \ [a, b] where H 2 is the upper half-plane {z ∈ C : Im (z) > 0} and, a and b are two complex numbers in H 2 (see Figure 18 (left)).
The method presented in Section 2 is not directly applicable to such domains Ω. So, we first map this domain to a domain G of the forms considered in Section 2. Since there is no exact conformal mapping from the domain Ω onto a doubly connected domain G bordered by smooth Jordan curves, we find such equivalent domain G using numerical methods. An iterative numerical method for computing such domain G has been presented in [NG]. We will omit the details here about the iterative method and refer the reader to [NG]. Using this iterative method, for a give half-plane with a segment domain Ω, a conformally equivalent domain G bordered by smooth Jordan curves can be obtained as in Figure 18 (right). For the new domain G, we use the MATLAB function annq with n = 2 11 to compute the capacity of the given domain Ω.
We tested our methods for several values of s and r. First, we chose the vertical segment F = [si, (1 + s)i], i.e., a = si and b = (1 + s)i, for 0.05 ≤ s ≤ 6. The relative errors of the computed values of the capacity for this case are presented in Figure 19. We see from Figure 19 that the presented method gives accurate results with relative error around 10 −14 . The approximate values of the capacity, the exact values of the capacity, and the total CPU time for several values of s and r are presented in Table 7. We also compute the values of the capacity for the vertical segment F = [(3−s)i, (3+s)i] for 0.05 ≤ s ≤ 2.95 and for the horizontal segment F = [−s + 3i, s + 3i] for 0.05 ≤ s ≤ 3. Both segments pass through the point 3i and have the length 2s. The results are presented  in Figure 20. As we can see from Figure 20, the capacity increases as the length of the segment increases. For vertical segment, the capacity increases more rapidly when the segment becomes close to the real line. Finally, for a given point z 1 in H 2 , we define the function u(x, y) by for −3 < x < 3 and 0 < y < 3 such that x + iy = z 1 . We plot the contour lines for the function u(x, y) corresponding to several levels. The contour lines are shown in Figures

A strip with a slit
In this section, we consider the doubly connected domain Ω = S \ F where S is the infinite strip |Im z| < π 2 , a, b ∈ S , and F is the segment [a, b] (see Figure 23(left)). Similar to the domains considered in Sections 4 and 6, the method presented in Section 2 is not directly applicable the domain Ω = S \ F describe above. Thus, we first map this domain onto an equivalent domain G of the forms considered in Section 2, i.e., a domain bordered by smooth Jordan curves ((see Figure 23(right)).
For the two domains consider in the Sections 4 and 6, we used the iterative method presented in [NG] to find equivalent domains bordered by smooth Jordan curves. Similar iterative methods have been presented in [N4] for other canonical slit domains. However, the canonical domain Ω = S \ F obtained by removing a segment F from the infinite strip S has not been considered in [NG,N4]. So, for the domain Ω, we shall present in the next subsection an iterative method for computing an equivalent doubly connected domain G bordered by smooth Jordan curves. The iterative method is similar to the iterative methods presented in [AST,NG,N4]. 7.1. A preimage domain for a strip with a slit domain. Let G be the doubly connected domain in the interior of the unit circle Γ 1 and the exterior of a smooth Jordan curve Γ 2 . Let also Ω be the canonical domain consisting of the strip |Im z| < π/2 with a slit L along straight line such that the angle between the line and the positive real axis is θ 2 ; see Figure 23 (left). An efficient numerical method for computing the conformal map z = Φ(ζ) from G onto the domain domain Ω such that has been presented in [NF]. The center and the length of the slit L depend uniquely on the given domain G. Assume that Γ = ∂G is parametrized by a function η(t) as in (2.5). Assume also θ 1 = 0, θ 2 is the angle between the slit L and the positive real axis, the function A be defined by (2.6); and the operators N and M are as in Subsection 2.3. Then, the following theorem from [NF] provides us with a method for computing the conformal mapping from the domain G onto the domain Ω.
Theorem 7.3. Let θ 1 = 0, let θ 2 be the angle between the slit L and the positive real axis, let the function A be defined by (2.6), and let the function γ be defined by where the function Ψ is defined by If ρ is the unique solution of the boundary integral equation (2.9) and the piecewise constant function h is given by (2.10), then the conformal mapping Φ from G onto Ω is given by where f is the analytic function in G with the boundary values (2.11).
In Theorem 7.3, the domain G is assumed to be given and the canonical domain Ω is to be found. The integral equation (2.9) is then used to find the domain Ω and the conformal map z = Φ(ζ) from G onto Ω. In our application, however, we assume that the domain Ω (strip with slit) is given and we want to find the domain G. For such case, a direct application of Theorem 7.3 is not possible. So, similar to the iterative methods presented in [AST,NG,N4], we will describe in this subsection an iterative method to compute the (unknown) preimage domain G and the conformal map z = Φ(ζ) from the domain G onto the (known) strip with a slit domain Ω such that Φ satisfies the normalization conditions (7.2).
For a given doubly connected domain Ω consisting of the strip |Im z| < π/2 with a rectilinear slit L, let |L| denote the length of the slit L, ζ L denotes its center, and θ 2 denotes the angle between the positive real axis and the slit L. For k = 0, 1, 2, 3, . . ., where k denotes the iteration number, letΩ k be the doubly connected connected domain in the strip |Im ξ| < π/2 and in the exterior of the ellipseL k parametrized bŷ η k (t) =ẑ k + 0.5a k e iθ 2 (cos t − ir sin t), t ∈ J 2 , where is the ratio of the lengths of the major to the minor axes of the ellipse, 0 < r ≤ 1. The parametersẑ k and a k of the ellipses will be computed by the iterative method described below. Note that, the domainΩ k is obtained from Ω by replacing the slit L by a thin ellipsê L k whose major axis is on the slit L (see Figure 23 (center)). Then the preimage domain G k is assumed to be the image of the doubly connected domainΩ k under the conformal mapping ζ = Ψ −1 (ξ) where Ψ −1 is the inverse of the function Ψ in (7.5), i.e., (see Figure 23 (right)). Hence, the preimage domain G k is the bounded doubly connected domain interior to the unit circle parametrized by η k 1 (t) = e it , t ∈ J 1 , and exterior to the quasi-ellipse Γ 2 parametrized by η k 2 (t) = Ψ −1 η k (t) , t ∈ J 2 . The parametersẑ k and a k will be computed as follows: Initialization: Setẑ 0 = ζ L , a 0 = (1 − 0.5r)|L|, j = 1, 2, . . . , m.

Iterations:
For k = 1, 2, 3, . . ., • The preimage domain G k−1 is mapped by the method presented in Theorem 7.3 onto the canonical domain Ω k (the strip |Im ζ| < π/2 with a rectilinear slits L k making angles θ 2 with the positive real axis). • Let ζ k L be the center of the slit L k and |L k | be its length. The parameters of the preimage domain G k are updated througĥ where Max is the maximum number of allowed iterations and ε is a given tolerance. In our numerical experiments, we choose Max = 100 and ε = 10 −14 . The iterative method produces a sequence of doubly connected domains G 0 , G 1 , G 2 , G 3 , . . . which converges to the required preimage domain G bordered by smooth Jordan curves. The iterative method provides us also with a conformal map z = Φ(ζ) from G onto the given domain Ω.
Finally, similar iterative methods have been tested in [NG,N4] for other canonical domains. The numerical examples presented in [NG,N4,LSN] show that the iterative method converges after only few number iterations even for domain with high connectivity. However, so far, no proof of convergence is yet available. 7.7. Computing the capacity. For a given domain Ω = S \ F , in this example, we use the above iterative method with n = 2 11 to find an equivalent domain G bordered by smooth Jordan curves. Then we use the MATLAB function annq to compute the capacity of G which is equal to the capacity of the given domain Ω. For the segment F = [0, s] where s > 0 is a real number, we can find the exact value the capacity of Ω. It is clear that the mapping w = Φ(z) = e z − 1 e z + 1 = th(z/2) maps conformally the strip − π 2 < Im z < π 2 onto the unit disk |w| < 1 and maps the segment F = [0, s] onto the segment [0, tanh(s/2)] on the real line interior to the unit circle. Then capacity of the domain obtained by removing the real segment [0, tanh(s/2)] from the unit disk is known and is given by (see [LV], [V2,Thm 8.6(1)]) (7.8) 2π µ(tanh s 2 ) which is also the capacity of the domain Ω by conformal invariance of the capacity. For such a domain, we use our method to compute the capacity for several values of s. The approximate values of the capacity, the exact values of the capacity, and the total CPU time are presented in Table 8.
The presented method is used also to compute approximate values of the capacity for other segments. The computed values of the capacity for the vertical segment F = [(s − 0.5)i, (s + 0.5)i] for −1.05 ≤ s ≤ 1.05 vs s are shown in Figure 24 (left). Figure 24 (left) shows also the computed values of the capacity for the horizontal segment F = [−0.5 + si, 0.5 + si] for −1.55 ≤ s ≤ 1.55. As we can see from Figure 24 (left), the capacity increases as the segment moves vertically toward the lines Im z = −π/2 and Im z = π/2. Figure 24 (right) shows the computed values of the capacity for the vertical segment F = [s − 0.5i, s + 0.5i] and the horizontal segment F = [s − 0.5, s + 0.5] for −4 ≤ s ≤ 4. As we can see from Figure 24 (right), the capacity does not change as the segment moves horizontally in parallel to the lines Im z = −π/2 and Im z = π/2. We also compute the values of the capacity for the vertical segment F = [−si, si] and the horizontal segment F = [−s, s] for 0.05 ≤ s ≤ 1.5. Both segments passes through the origin and has the length 2s. The results are presented in Figure 25. As we can see from Figure 20, the capacity increases as the length of the segment increase. For vertical segment, the capacity increases more rapidly when the segment becomes close to the lines Im z = −π/2 and Im z = π/2. Finally, for a given point z 1 in S, we define the function u(x, y) by for −3 < x < 3 and −π/2 < y < π/2 such that x + iy = z 1 . We plot the contour lines for the function u(x, y) corresponding to the several levels. The contour lines are shown in Figures 26 for z 1 = 0 and in Figures 27 for z 1 = i. Table 9 presents the values of the capacity cap(S\[z 1 , z 2 ]) for several values of z 1 and z 2 .   where 0 < d < 0.5 (see the Figure 28). We use the MATLAB function annq presented in Subsection 2.23 with n = 2 15 to compute the capacity of G for several values of d.
When d = 0, the two rectangles reduced to the two slits [i/2, 1 + i/2] and [−i/2, 1 − i/2]. For these two slits, we can use the numerical method presented in Section 4 to compute the capacity of the domain in the exterior to these two slits. The obtained results are presented at the bottom of Table 10.
Let R 1 be the unbounded doubly connected domain in the exterior to the two slits  domain R 2 equals to the capacity of R 3 which can be expressed by [V2, 5.60 (1)] Here µ is the function defined in (3.7).
By [BF,119.03], the domain R 2 can be mapped conformally also onto the unbounded doubly connected domain R 4 in the exterior of the two slits [−t/2 − ib/2, −t/2 + ib/2] and [t/2 − ib/2, t/2 + ib/2] with where the functions E(k), K(k) are defined in (3.6) and (3.5), resp., and Hence cap(R 4 ) = cap(R 2 ). Further, it is clear that the domain R 1 can be conformally mapped by the functionΨ(z) = i(z − 1/2) onto the domain R 4 if we choose k such that t = b = 1. Thus, the exact value of the capacity of R 1 is given by where k satisfies the equations The equations (8.4) are solved using Mathematica for k and the value of the capacity of R 1 computed through (8.3) is 2.1157789709245134. This value agrees with the value presented at the bottom of Table 10 with relative error 1.5 × 10 −9 . 8.5. A vertical rectangle in the upper half-plane. Consider the doubly connected domain G in the exterior of the rectangular closed set 1,2] in the upper half-plane where 0 < d < 0.5 (see the Figure 29 for d = 0.1).
The the auxiliary map (8.6) w = Ψ(z) = iz + 1 z + i in the upper half-plane where 0 < d < 0.5 (see the Figure 30). By symmetry, the capacity for this domain is 2 times the capacity for the two rectangles case considered in Subsection 8.1. As in the previous example, the the auxiliary map w = Ψ(z) in (8.6) is used to transform the domain G onto a domainĜ interior to the unit disk and exterior to a piecewise smooth Jordan curve (see Figure 30). Then G andĜ have the same capacities. We use the function annq with n = 2 15 to compute the capacity ofĜ for several values of d.
When d = 0, the rectangle reduced to the slit [0.5i, 1 + 0.5i]. By symmetry, the capacity for the half-plane with the horizontal slit [0.5i, 1+0.5i] is 2 times the capacity for the cases of the domain exterior to the two horizontal slits [i/2, 1 + i/2] and [−i/2, 1 − i/2] considered in Subsection 8.1. Thus, according to the exact capacity presented in Subsection 8.1, the exact capacity for the upper half-plane with the horizontal slit [0.5i, 1+0.5i] is 4.23155794184903.
For numerical computing of the capacity of the upper half-plane with the slit [0.5i, 1+0.i], we use the method presented in Section 6. The obtained result is presented at the bottom  Finally, the third column in Table 12 shows halves of the computed values of the capacity for the domain presented in this section. The values presented in the third column agrees with the results presented in Table 10 for two rectangle case. 8.8. A rectangle in a strip. Consider the doubly connected domain G in the exterior of the rectangular closed set [0, 1] × [−d, d] in the infinite strip |Im z| < π 2 where 0 < d < 0.5 (see the Figure 31). The the auxiliary map w = Ψ(z) = tanh z 2 is used to transform the domain G onto a domainĜ interior to the unit disk and exterior to the piecewise smooth Jordan curve L which is the image of the rectangle under the Figure 31. The domain G in the exterior of a rectangular closed set in the infinite strip |Im z| < π 2 (left) and its imageĜ under the auxiliary map Φ (right) for d = 0.1. is computed using the numerical method presented in Section 7. The obtained numerical results are presented in Table 13. For the domain exterior to slit [0, 1] in the strip |Im z| < π 2 , the exact value of the capacity can be computed from (7.8) and is equal to 2π/µ(tanh(1/2)) = 2.99266869365819. This exact value agrees with the result presented at the bottom of Table 11 with relative error 2.8 × 10 −13 . 9. The hyperbolic capacity and the elliptic capacity Let E be a compact and connected set (not a single point) in the unit disk D. In this section, we use the MATLAB function annq in Subsection 2.23 to compute the hyperbolic capacity and the elliptic capacity of the set E. Both the hyperbolic capacity and the elliptic capacity are invariants under conformal mappings 9.1. The hyperbolic capacity The hyperbolic capacity of E, caph(E), is defined by [Va,p. 19] (9.2) caph(E) = lim n→∞ max z 1 ,...,zn∈E 1≤k<j≤n .
For the hyperbolic capacity, we assume G is the bounded doubly connected domain defined by G = D\E such that G can be mapped conformally onto an annulus q < |w| < 1. Hence the hyperbolic capacity caph(E) is given by [DK] (9.3) caph(E) = q.
The constant q can be computed by the function annq.
9.4. The elliptic capacity For the compact and connected set E, we define the antipodal set E * = {−1/a : a ∈ E}. Since we assume E ⊂ D, we have E ∩ E * = ∅ (in this case, the set E is called "elliptically schlicht" [DK]). The elliptic capacity of E, cape(E), is defined by [DK] (9.5) cape(E) = lim n→∞ max z 1 ,...,zn∈E 1≤k<j≤n .
To compute the elliptic capacity, we assume G is the doubly connected domain between E and E * such that G can be mapped conformally onto an annulus r < |w| < 1/r. Then the elliptic capacity is given by [DK] cape(E) = r.
Here, the domain G could be bounded or unbounded. We shall use the method presented in Section 2 to map the domain G onto an annulus q < |w| < 1 which is conformally equivalent to the annulus r < |w| < 1/r with r = √ q. Thus, we have (9.6) cape(E) = √ q.
We compute q using the function annq. Finally, as our interest in this paper is only in closed and connected subset E of the unit disk D and comparing numerically between the values of cape(E) and caph(E), it is worth mentioning that Duren and Kühnau [DK] have proved that with equality if and only if E = −E. This inequality is verified numerically in the following numerical examples. 9.7. A disk. As our first example, we compute the hyperbolic capacity and the elliptic capacity of the disk E = {z ∈ C : |z| ≤ r}, 0 < r < 1. For this set E, both capacities are equal where [Ki, DK] caph(E) = cape(E) = r. For computing caph(E), we use the function annq with α = (1 + r)/2 and z 2 = 0 to compute the value of q for the conformal map of the double connected domain G = D\E (see Figure 32 (left)) onto the annulus q < |w| < 1 and hence caph(E) = q. For cape(E), the domain G between E and E * is the bounded doubly connected domain r < |z| < 1/r (see Figure 32 (right)). We use the MATLAB function annq with α = 1 and z 2 = 0 to compute the value of q for the conformal map of this domain G onto the annulus q < |w| < 1 and hence cape(E) = √ q. For both cases, we use n = 2 12 and 0.02 ≤ r ≤ 0.98. The relative error in the obtained results for caph(E) and cape(E) are shown in Figure 33.
9.8. A square. For the second example, we assume E is the closed set [−r, r] × [−r, r], 0 < r < 1/ √ 2. For computing caph(E), the domain G is the bounded doubly connected domain in the interior of the unit circle and in the exterior of the square (see  (left)). We use the function annq with α = (1 + r)/2 and z 2 = 0 to compute q and then caph(E) = q. For cape(E), the domain G is the bounded doubly connected domain between E and E * (see Figure 34 (right)). Hence, cape(E) = √ q where q is computed using the function annq with α = (r + 1/r)/2 and z 2 = 0. For both cases, we use n = 2 13 for 0.02 ≤ r ≤ 0.69. The obtained results are shown in Figure 35. This set E is symmetric where E = −E, and hence caph(E) = cape(E).
9.9. Amoeba-shaped boundary. For the third example, we compute caph(E) and cape(E) of E where E is the closed region bordered by the amoeba-shaped boundary L with the parametrization η(t) = 0.1 + 0.6i + 0.2 e cos t cos 2 2t + e sin t sin 2 2t e −it , 0 ≤ t ≤ 2π.
For the hyperbolic capacity caph(E), the domain G is the bounded doubly connected domain in the interior of the unit circle and in the exterior of the curve L (see Figure 36 (left)). Then caph(E) = q where q is computed using the function annq with α = 0 and z 2 = 0.25 + 0.5i. To compute cape(E), the domain G is the unbounded doubly connected domain in the exterior of E and E * (see Figure 36 (right)). We use the function annq with  The exact value of the hyperbolic capacity of E is [Ki] caph(E) = e −µ(r) .
For computing caph(E) numerically, the domain G is the domain exterior of the segment [0, r] and interior to the unit circle |z| = 1. We use the the mapping function ζ = (Ψ −1 2 • Ψ 1 )(z) to conformally map the domain G onto the domainĜ exterior to the the unit circle |ζ| = 1 and interior to a smooth Jordan curve which is the image of the unit circle |z| = 1 under the mapping function Ψ −1 2 • Ψ 1 . Then, we use the function annq with α = Φ −1 2 (Φ 1 (0.5i)) and z 2 = 0 to compute the value of q for the conformal map of this domainĜ onto the annulus q < |w| < 1 and hence, by (9.3), caph(E) = q.
To use the the presented method to compute the elliptic capacity cape(E) numerically, we use the mapping function ζ = Ψ 4 (z) = z τ , to map the the bounded domain G 2 onto the bounded domain G 3 exterior to the segment [0, 1] and interior to the circle |ζ| = 1/τ . Finally, the function ξ = Ψ −1 2 (ζ) maps the domain G 3 onto the bounded domainĜ exterior to the unit circle and interior to as smooth Jordan curve which is the image of the circle |ζ| = 1/τ under the mapping function Ψ −1 2 . Thus, the bounded doubly connected domainĜ is the image of the original domain G under the mapping Ψ −1 2 • Ψ 4 • Ψ 3 • Ψ −1 2 • Ψ 1 . Then, we use the function annq with α = Φ −1 2 (Φ 4 (0.5i)) and z 2 = 0 to compute the value of q for the conformal map from this domainĜ onto the annulus q < |w| < 1 and hence, by (9.6), cape(E) = √ q. For both cases, the relative error in the obtained results for 0.02 ≤ r ≤ 0.98 is shown in Figure 37 (right). Figure 37 (left) shows also the computed capacities. For both cases, we use n = 2 12 .

Concluding Remarks
Conformal invariants are important tools for complex analysis with many applications. However, these invariants can be expressed explicitly only in very few special cases. Thus, numerical methods are required to compute these invariants. A numerical method for computing some conformal invariants is presented in this paper. The method can be used for domains with different types of boundaries including domains with smooth or piecewise smooth boundaries. The performance and the accuracy of the presented method is compared to analytic solutions or to previous results whenever analytic solutions or previous results are available. Further, a MATLAB implementation of the proposed method is given in the MATLAB function annq in Subsection 2.23. This MATLAB function was used in almost all examples in this paper to compute the conformal capacity, the hyperbolic capacity and the elliptic capacity. For some examples, an auxiliary procedure is required before using the function annq. All the computer codes of our computations are available in the internet link https://github.com/mmsnasser/cci.