Conformal Invariants in Simply Connected Domains

This paper studies the numerical computation of several conformal invariants of simply connected domains in the complex plane including, the hyperbolic distance, the reduced modulus, the harmonic measure, and the modulus of a quadrilateral. The used method is based on the boundary integral equation with the generalized Neumann kernel. Several numerical examples are presented. The performance and accuracy of the presented method is validated by considering several model problems with known analytic solutions.

originally proven for functions defined in the unit disk to the case when the domain of definition is a simply connected domain. Therefore for the convenient analysis of distances and other metric characteristics of sets it is natural to use conformally invariant distances and set functions. This works fine in the cases when the Riemann mapping function is known explicitly, such as in the cases described in [12]. Polygons form a large class of plane domains for which the Riemann mapping function can be given in terms of the Schwarz-Christoffel formula which is semi-explicit, although it involves unknown accessory parameters. A numerical implementation of conformal mapping methods based on the Schwarz-Christoffel formula is documented in [6] and the Schwarz-Christoffel Toolbox [5] is widely used to solve mapping problems in concrete applications. The so called crowding phenomenon, an inherent computational challenge in these mapping problems, is described in [6, pp. 20-21], [21, pp. 75-77]. This phenomenon can be observed already in numerical conformal mapping of rectangles onto a half-plane when the quotient m > 1 of the side lengths is large enough. The critical value of m depends on the computer floating point arithmetic and for double precision arithmetic the critical value lies in the range [10,20] [6, pp. 20-21], [21, pp. 75-77].
We apply here the boundary integral equation method as developed in [14,16,17,19,25] to compute numerically conformal invariants such as the hyperbolic metric, harmonic measure, reduced modulus, and the modulus of a quadrilateral [7,8,21,23]. The cases considered here include, in particular, classes of domains to which the earlier methods do not seem to apply. Our methods are described in the respective sections of the paper, they are implemented in MATLAB, and the results are summarized by tabular data and graphics. We also give experimental error estimate in some simple cases. We include some code snippets within the text to indicate implementation details. All the computer codes of our computations are available in the internet link https://github.com/mmsnasser/ci-simply. Section 2 reviews the boundary integral equation method for computing the conformal mapping from bounded and unbounded simply connected domains onto circular domains. This method will be applied in the remaining sections, sometimes together with auxiliary procedures. In Sect. 3 we use our method to compute the hyperbolic distance for several examples. Section 4 deals with the reduced modulus for bounded and unbounded simply connected domains. Section 5 deals with the harmonic measure for a simply connected polygonal domains. In Sect. 6, we present an iterative method for numerical computation of the conformal mapping from a quadrilateral onto a rectangle. We also present an analytic example to illustrate the effect of crowding phenomenon on the accuracy of such mapping.

Conformal Mappings of Simply Connected Domains
In this section, we review a numerical method for computing the conformal mapping from bounded and unbounded simply connected domains onto the unit disk and the exterior unit disk, respectively. The method is based on the boundary integral equation with the generalized Neumann kernel [16][17][18]25].

The Generalized Neumann Kernel
Let G be a bounded or an unbounded simply connected domain bordered by a closed smooth Jordan curve Γ = ∂G. The orientation of the boundary Γ is counterclockwise when G is bounded and clockwise when G is unbounded. We parametrize Γ by a 2πperiodic complex function η(t), t ∈ [0, 2π ]. We assume that η(t) is twice continuously differentiable with η (t) = 0 (the presented method can be applied also if the curve Γ has a finite number of corner points but no cusps [19]). We denote by H the space of all Hölder continuous real functions on the boundary Γ .
Let A be the complex function [17] The kernel N (s, t) is continuous [25]. Hence, the integral operator N defined on H by is compact. The integral equation with the generalized Neumann kernel involves also the following kernel which is singular and its singular part involves the cotangent function [25]. The integral operator M defined on H by is singular, but it is bounded on H [25].

Bounded Simply Connected Domain
Let w = Φ(z) be the unique conformal map of the bounded simply connected domain G onto the unit disk |w| < 1 such that Then, the mapping function Φ with normalization (4) can be written for z ∈ G ∪ Γ as [18, § 3] where the function f (z) is analytic in G with the boundary values ρ is the unique solution of the integral equation and c = e −h where the constant h is given by Notice that Φ (α) = c = e −h > 0. Instead of the normalization (4), we can assume that the mapping function Φ satisfies the normalization In this case, the function Φ maps the domain G onto the disk |w| < R with a positive constant R. The constant R is uniquely determined by G and the point α and is called the conformal radius of G with respect to α and is denoted by R(G, α). For this case, in view of (5), we can write the mapping function Φ for z ∈ G ∪ Γ as where the function f is as in (5), i.e., we divide the right-hand side of (5) by c = e −h . Hence, the boundary values of the mapping function Φ satisfy |Φ(η(t))| = 1/c = e h which implies where the constant h is as in (7).

Unbounded Simply Connected Domain
For an unbounded simply connected domain G ⊂ C with ∞ ∈ G, there exists a unique conformal map w = Φ(z) from G onto the exterior of the unit disk |w| > 1 such that Then, the mapping function Φ with normalization (4) can be written for z ∈ G ∪ Γ as [18, § 3] where β is an auxiliary point in G c = C\G and f (z) is an analytic function in G with f (∞) = 0. The boundary values of the function f are given by where A(t) = 1, the function γ is defined by γ (t) = − log |η(t) − β|, ρ is the unique solution of the integral equation (6), and c = e −h where the constant h is given by (7). Notice that Φ (∞) = c = e −h > 0. Similar to the bounded case, the normalization (11) can be replaced by For the new normalization (13), the mapping function Φ maps the unbounded domain G onto the exterior disk |w| > R with a positive constant R. The constant R is uniquely determined by G and is called the conformal radius of G with respect to ∞ and is denoted by R(G, ∞). If Φ satisfies the normalization (13), then it can be written as and its boundary values satisfy |Φ(η(t))| = 1/c = e h where the function f is as in (5) and h is given by (7). Hence,

Numerical Solution of the Boundary Integral Equation
The integral equation (6) is solved using the MATLAB function fbie in [17]. The function fbie is based on using the MATLAB function zfmm2dpart in the toolbox In the numerical experiments in this paper, the parameters in fbie are chosen as following: iprec = 5 (the tolerances of the FMM is 0.5 × 10 −15 ), gmrestol = 0.5 × 10 −14 (the tolerances of the GMRES), restart = [ ] (the GMRES is used without restart), and maxit = 100 (the maximum number of iterations for GMRES). Finally, the values of the auxiliary points α in (5), (9) and β in (12), (14) have no effects on the values of the conformal mapping Φ as long as these points are sufficiently far away from the boundary Γ .

Hyperbolic Distance
For x, y ∈ D the hyperbolic distance ρ D (x, y) is defined by [4,8,11] sinh The main property of the hyperbolic distance is the invariance under Möbius transformations of D onto D defined by where a ∈ D is fixed. In the metric space (D, ρ D ) one can build a non-euclidean geometry, where the parallel axiom does not hold. In this geometry, usually called the hyperbolic geometry of the Poincare disk, lines are circular arcs perpendicular to the boundary ∂D . This geometry is fundamentally different from the Euclidean geometry, but many results of Euclidean geometry have counterparts in the hyperbolic geometry [4].
Let G be a Jordan domain in the plane. One can define the hyperbolic metric on G in terms of the conformal Riemann mapping function Φ : G → D = Φ(G) as follows: This definition yields a well-defined metric, independent of the conformal mapping Φ [4,8,11]. In hyperbolic geometry the boundary ∂G has the same role as the point of {∞} in Euclidean geometry: both are like "horizon", we cannot see beyond the horizon. It turns out that the hyperbolic geometry is more useful than the Euclidean geometry when studying the inner geometry of domains in geometric function theory.

Computing the Hyperbolic Distance for Simply Connected Domains
Let G ⊂ C be a bounded simply connected domain, let α ∈ G, and let w = Φ(z) be the unique conformal map of G onto the unit disk |w| < 1 with the normalization (4). Then for any two points z 1 and z 2 in G, we can define the hyperbolic metric ρ G in the usual way, (16) A MATLAB function for computing the hyperbolic metric ρ G (z 1 , z 2 ) for any two points z 1 and z 2 in the bounded simply connected domain G is given below.  *(1 -abs ( Phio ) .^2) ) ) ; end In the remaining part of this section, we use the MATLAB function hypdist to compute the hyperbolic metric ρ G (z 1 , z 2 ) for several examples. In these examples, for a given point z 1 in G, we define the function u(x, y) for all x and y such that Then we use the MATLAB function hypdist to compute the values of the function u(x, y) in the domain G and plot the contour lines for the function u(x, y) corresponding to the several levels.

L-Shaped Polygon
As our first example, we consider the simply connected domain G in the interior of the L-shaped polygon with the vertices 6 + i, 1 + i, 1 + 4i, −1 + 4i, −1 − i, and 6 − i. The contour lines of the function u(x, y) obtained with n = 6 × 2 9 discretization points on the boundary of the L-shaped polygon are shown in Fig. 1 (left) for z 1 = 2i and in Fig. 1 (right) for z 1 = 2. Table 1 presents the values of the hyperbolic metric ρ G (z 1 , z 2 ) for z 1 = 2i (left), z 1 = 2 (right), and for several values of z 2 . Table 1 The values of the hyperbolic metric ρ G (z 1 , z 2 )

Fig. 2
The contour lines of the function u(x, y) for the amoeba-shaped boundary for z 1 = 2 (left) and

Amoeba-Shaped Boundary
In the second example, we consider the simply connected domain G in the interior of the curve Γ (amoeba-shaped boundary [3]) with the parametrization η(t) = e cos t cos 2 2t + e sin t sin 2 2t e it , 0 ≤ t ≤ 2π.
The contour lines for the function u(x, y) computed with n = 2 12 are shown in Fig. 2 (left) for z 1 = 2 and in Fig. 2 (right) for z 1 = −1 + i.

Circular Arc Quadrilateral
For the third example, we consider the simply connected domain G in the interior of the circular arc quadrilateral consisting of the right-half of the circle with center 1 and radius 1, the upper-half of the circle with center i and radius 1, the left-half of the circle with center −1 and radius 1, and the lower-half of the circle with center −i and radius 1. The contour lines for the function u(x, y) computed with n = 2 12 are shown in Fig. 3 (left) for z 1 = 1.5 and in Fig. 3 (right) for z 1 = 0.

Circular Arc Polygon
In the fourth example, we consider the simply connected domain G in the interior of the circular arc polygon with seven sides. The contour lines for the function u(x, y)

Reduced Modulus
The reduced modulus for simply connected domains is defined in terms of the conformal radius of simply connected domains introduced in Sect. 2.
For an unbounded simply connected domain G ⊂ C with ∞ ∈ G, the reduced modulus of the domain G with respect to ∞ is defined by [23, p. 17] where R = R(G, ∞) is the conformal radius of G with respect to ∞.

Computing the Reduced Modulus of Simply Connected Domains
As was explained in Sect. 2, the conformal radius of simply connected domains can be computed using the integral equation with the generalized Neumann kernel. For bounded simply connected domains, it follows from (10) that the reduced modulus of the domain G with respect to the point α is given by where the constant h is given by (7). For unbounded simply connected domains, it follows from (15) that the reduced modulus of the domain G with respect to the point ∞ is given by The above method for computing the conformal radius and the reduced modulus of bounded and unbounded simply connected domains can be implemented in MATLAB as in the following function.

Domain Exterior to an Ellipse
As our first example, we consider the simply connected domain G r in the exterior of the ellipse For r = 0, the ellipse reduces to the segment [−1, 1] and for r = 1 to the unit circle.
We can easily show that the function 1 w maps the domain exterior to the circle |w| = (1 + r )/2 onto the domain exterior of the ellipse. Hence, the inverse mapping maps the domain exterior to the ellipse onto the domain |w| > (1 + r )/2, where the branch of the square root is chosen such that We use the MATLAB function confrad to compute the reduced modulus m(G r , ∞) with n = 2 12 for 0.005 ≤ r ≤ 1. The obtained results are shown in Fig. 5.

Domain Interior to an Ellipse
For the second example, we consider the simply connected domain G r in the interior of the ellipse Let w = Φ(z) be the unique conformal mapping from the interior of the ellipse onto the interior of the unit circle with the normalization Φ(0) = 0 and Φ (0) > 0. The exact form of the inverse conformal mapping z = Φ −1 (w) is given in [10]. In particular, it was shown in [10]  where μ −1 is the inverse of the Grötzsch modulus function, see (28) below. Hence, Φ (0) = 2 √ s K (s)/π . Thus, the mapping functionΦ defined by is the unique conformal mapping from the interior of the ellipse onto the with the normalizationΦ(0) = 0 andΦ (0) = 1. Thus, R(G r , 0) = π/(2 √ s K (s)) and hence We use the MATLAB function confrad to compute the reduced modulus m(G r , 0) with n = 2 12 for 0.2 ≤ r ≤ 20. The obtained results are shown in Fig. 6.

Slitted Unit Disk
In the third example, we consider three types of slitted unit disks. In each case the exact reduced moduli are given in [23, p. 33].
1) G 1 = D\(−1, 0] where D is the unit disk. The exact value of the reduced modulus of G 1 with respect to r ∈ (0, 1) is given by [23] To use the integral equation to compute m(G 1 , r ), we first use the auxiliary map where the branch of the square root is chosen on the negative real line, to open up the slit and map the region G 1 onto a regionĜ 1 bordered by piecewise smooth Jordan curve where Φ 1 (r ) = 2r and Φ 1 (r ) = 1.  Fig. 7 (left).
2) G 2 = D\[r , 1) for 0 < r < 1. The exact value of the reduced modulus of G 2 with respect to the origin is given by [23] To compute m(G 2 , r ), we first use the auxiliary map where the branch of the square root is chosen on the positive real line, to map the region G 2 onto a regionĜ 2 bordered by piecewise smooth Jordan curve where Φ 1 (0) = −2r and Φ 1 (0) = 1. Then, m(G 2 , 0) = m(Ĝ 2 , −2r ). We use the MATLAB function confrad to compute m(Ĝ 2 , −2r ) with n = 2 12 for 0.01 ≤ r ≤ 0.99. The obtained results are shown in Fig. 7 (center).
The exact value of the reduced modulus of G 3 with respect to r ∈ (a, 1) is given by [23] To compute m(G 3 , r ), we first use the auxiliary map where the branch of the square root is chosen on the negative real line, to map the region G 3 onto a regionĜ 3 bordered by a piecewise smooth Jordan curve where

Polygon
For the fourth example, we consider the simply connected domain G in the interior of a polygon with vertices where ≥ 3 (see Fig. 8 (left) for = 8). We assume that the vertices of the polygon are given by v k = e 2kπ i/ , k = 0, 1, 2, . . . , − 1.
In this example, the exact value of the reduced modulus is unknown. We use the MATLAB function confrad to compute the reduced modulus m(G , 0) with n = × 2 9 for = 3, 4, . . . , 40. The obtained results are shown in Fig. 8. It is clear from this figure that m(G , 0) < 0 which means that R(G , 0) < 1 for the above values of . In other words, the conformal mapping Φ with the normalization (8) maps the domain G onto a disk interior to the unit disk.

Harmonic Measure
Let G be a Jordan domain in C and Γ be its boundary. Let also L be a boundary arc on Γ such that L = ∅ and Γ \L = ∅. The harmonic measure of L with respect o G is the C 2 (G) function u : G → (0, 1) satisfying the Laplace equation in G and u(z) → 1 when z → L and u(z) → 0 when z → Γ \ L . The harmonic measure is one of the key notions of potential theory and it has numerous applications to geometric function theory [8]. The harmonic measure of L with respect to G will be denoted by ω(z, L) (see e.g., [

Harmonic Measure for the Unit Disk
Assume that G is the unit disk |z| < 1, Γ is the unit circle |z| = 1, and L is the right half of the unit circle. It is clear that the Möbius transformation where the branch with log 1 = 0 is chosen.
To compute ω(z, L), we discretize the parametrization η(t), 0 ≤ t ≤ 2π , of the polygon Γ on each segment [z k , z k+1 ] by n s graded points on [2(k − 1)π/m, 2kπ/m]. Thus, the whole polygon Γ is discretized by n = mn s point t i , i = 1, 2, . . . , n in [0, 2π ] such that z k = η(t 1+(k−1)n s ) for k = 1, 2, . . . , m. Then we use the method presented in Sect. 3 to compute the conformal mapping ζ = Φ(z) from the interior of Γ onto the unit disk |ζ | < 1. The mapping function Φ maps the two points z k and z k+1 onto two points ζ 1 and ζ 3 , respectively, on the unit circle |ζ | = 1. The segment L is then mapped onto the arcL on the unit circle |ζ | = 1 from ζ 1 to ζ 3 . Let ζ 2  Fig. 9 The arc L between z 1 and z 3 and be the point on the middle ofL between ζ 1 and ζ 3 so that ζ 1 , ζ 2 and ζ 3 arranged in counterclockwise orientation [(see Fig. 9 (center)]. Then the Möbius transformation maps the unit disk |ζ | < 1 onto the unit disk |w| < 1 and maps the unit circle |ζ | = 1 onto the unit circle |w| = 1 such that the points ζ 1 , ζ 2 and ζ 3 are mapped onto the points −i, 1 and i, respectively. Thus, the mapping function Ψ maps the arcL on |ζ | = 1 onto the right half of the unit circle |w| = 1 [see Fig. 9 (right)]. Finally, the mapping function maps the domain G onto the disk |w| < 1 and maps the segment L on Γ onto the right half of the unit circle |w| = 1. Hence, by (21), the harmonic measure of L with respect to G is given by The above method for computing the harmonic measure of a segment L = [z k , z k+1 ] with respect to the polygon domain G can be implemented in MATLAB as in the function hm.m where the discretization of the parametrization of the polygon is computed using the MATLAB function polygonp.m (both functions can be downloaded from https://github.com/mmsnasser/ci-simply/tree/master/harmonic

Polygon with 5 Sides
As our first example, we consider the simply connected domain G in the interior of the polygon with 5 sides (the polygon shown in Fig. 10 and the vertices of the polygon are

Polygon with 13 Sides
For the second example, we consider the simply connected domain G in the interior of the polygon with 13 sides where the vertices of the polygon are 4, 4 + 2i, 2 + 4i, 4i, −1 + 3i, −2 + 3i, −3 + i, −3, −2 − 2i, −1 − 3i, −3i, 1 − 2i, and 3 − 2i. The MATLAB function hm with n s = 2 9 is used to compute the harmonic measure ω(z, L) for each side L of the polygon with respect to the polygon domain G. The level curves of the function ω(z, L) for the first 6 sides are shown in Fig. 11.

Quadrilateral Domains
Let w = Φ(z) be the conformal mapping from the interior of the unit circle D = {z ∈ C : |z| = 1} onto the interior of the rectangle R r = {w : 0 < Re w < 1, 0 < Im w < r } (22) such that If the domain R r is known (i.e., if r is known), then we can map R r onto the unit disk using the method described in Sect. 2. Let ζ = Ψ 1 (w) be the conformal mapping from R r onto the disk |ζ | < 1 such that where α is a given point in R r . The mapping function z = Ψ 1 (w) maps the points 0, 1, and 1 + ir on ∂ R r onto points ζ 1 , ζ 2 , and ζ 3 on the unit circle. Then the Möbius transform maps the unit disk |ζ | < 1 onto the unit disk |z| < 1 such that the points ζ 1 , ζ 2 and ζ 3 are mapped onto the points z 1 , z 2 and z 3 , respectively. Thus, the mapping function maps the domain R r onto the unit disk D such that the points 0, 1, 1 + ir are mapped onto the points z 1 , z 2 and z 3 , respectively. If z = Ψ (w) maps also the point ir onto

Iterative Method
In this section, for a given quadrilateral (D; z 1 , z 2 , z 3 , z 4 ), we present an iterative method for computing the unknown constant r and the mapping function z = Ψ (w) such that Ψ (0) = z 1 , Ψ (1) = z 2 , Ψ (1 + ir ) = z 3 , and Ψ (ir ) = z 4 . First we choose an initial value r 0 = 1, then we use the function Ψ to map R r 0 to a quadrilateral (D; z 1 , z 2 , z 3 , z 4,0 ) where z 4,0 is a point on the arc [z 3 , z 1 ] containing z 4 (see Fig. 12). The point z 4,0 could be on either side of z 4 on the arc [z 3 , z 1 ]. We add a correction Δ 0 to r 0 to get a new approximation r 1 . Then we map R r 1 to a quadrilateral (D; z 1 , z 2 , z 3 , z 4,1 ) using the function Ψ . The point z 4,1 is now close to the point z 4 . We continue with this iterative method to generate a sequence of approximation r 0 , r 1 , r 2 , . . . and the mapping function Ψ maps the rectangle R r k to a quadrilateral (D; z 1 , z 2 , z 3 , z 4,k ). We stop the iteration when the distance (on the unit circle) between the two points z 4 and z 4,k is small. Then, we consider r k as an approximation to r . Since, for each iteration k, the point z 4,k is on the arc [z 3 , z 1 ], we can always choose suitable corrections Δ k to ensure the convergence of the iterative method. In this paper, for k ≥ 1, we choose To accelerate the convergence of the iterative method, we introduce a factor δ k and we calculate r k using the formula where r 0 = 1, δ 0 = δ 1 = 1, and for k ≥ 2, (26) In other words, when the three points z 4,k−2 , z 4,k−1 and z 4,k are in the same side of z 4 , we double δ k−1 to increase the correction added to r k and so push z 4,k toward z 4 . However, when the three points z 4,k−2 , z 4,k−1 and z 4,k oscillate around z 4 , we bisect δ k−1 to reduce the correction added to r k . To avoid getting very long rectangle or very narrow rectangle during the iterations, we do not allow δ k Δ k to be more than 0.2r k or less than −0.2r k .

Algorithm
The above iterative method is summarized as follows. Initialization: Max is the maximum number of allowed iterations and ε is a given tolerance.
In our numerical experiments, we choose Max = 50 and ε = 0.5 × 10 −13 . The iterative method produces a sequence of numbers r 0 , r 1 , r 2 , r 3 , . . . which converges to the required constant r . The iterative method provides us also with a conformal map z = Ψ (w) from R r onto the given domain D. Then the required map Φ is given by The numerical examples presented in this section show that the iterative method converges for several examples. However, no proof of convergence is available so far.

Examples
We consider the computing of the conformal mapping from the quadrilateral domains (D; e iθ 1 , e iθ 2 , e iθ 3 , e iθ 4 ), onto rectangular domains for the following values of θ 1 , θ 2 , θ 3 , and θ 4 : The values of the modulus r = M(Q j ), j = 1, 2, 3, 4, the number of iterations required for convergence, and the total CPU time (sec) required for convergence are listed in Table 2. For the four domains, we use n = 2 11 . Orthogonal polar grids in the circular domains and their images under the conformal mapping are shown in Figs. 13, 14, 15 and 16. The points z 1 , z 2 , z 3 , z 4 on the unit circle and their images on the rectangle are shown as small colored squares such that a point z k and its image has the same color. For Q 3 , z 1 = e −0.5005π and z 2 = e −0.4995π which are very close to each other. Similarly, z 3 = e 0.4995π and z 4 = e 0.5005π are very close to each other. The length of the arcs between z 1 and z 2 and between z 3 and z 4 is 0.001π . Thus, we can not distinguish between z 1 and z 2 and between z 3 and z 4 in Fig. 15 (left). The small arc between z 1 and z 2 is mapped by the conformal mapping to the lower side of the rectangle. Similarly, the small arc between z 3 and z 4 is mapped by the conformal mapping to the upper side of the rectangle. Although these arcs are too small, the proposed iterative method converges after only 36 iterations. In Q 4 , the two points z 2 = e −0.0001π and z 3 = 1 are very close to each other where the length of the arcs between them is 0.0001π , and hence we can not distinguish between z 2 and z 3 in Fig. 16 (left). The small arc between z 2 and z 3 is mapped by the conformal mapping to the right side of the rectangle. The proposed iterative method converges after only 40 iterations.
For the three domains Q 2 , Q 3 , and Q 4 , the error per iteration is shown in Fig. 17. For Q 3 and Q 4 , we have points on the unit circle that are very close to each other. This explains why the number of iterations for Q 3 and Q 4 is larger than the number of iterations for Q 2 . For Q 1 , the method converges after only one iteration since the exact value of r is 1 which is the same as our initial value r 0 .  which can be mapped conformally onto the rectangular domain R r = {w : 0 < Re w < 1, 0 < Im w < r } such that the point 1 is mapped to 0, e iθ 1 is mapped to 1, e iθ 2 is mapped to 1 + r i, and e iθ 3 is mapped to r i. The quadrilateral domain can be mapped also by Möbius transform onto the upper half-plane such that the point 1 is mapped to −1, e iθ 1 is mapped to 0, e iθ 2 is mapped to a positive real number s, and e iθ 3 is mapped to ∞. See Fig. 18.
The exact value of the positive constant s can be obtained in terms of the values of θ 1 , θ 2 , and θ 3 . For distinct points z 1 , z 2 , z 3 , and z 4 in C, we define the absolute (cross) ratio by [4] This definition can be extended if z 4 = ∞ by taking the limit, i.e., Thus, for the four points 1, e iθ 1 , e iθ 2 , and e iθ 3 on the unit circle, we have Let Ω 1 be the family of curves lying in D and joining the arcs (1, e iθ 1 ) and (e iθ 2 , e iθ 3 ) (see Fig. 18). Similarly, let Ω 2 be the family of curves lying in R r and joining the segments (1, 1 + ir ) and (ir , 0), and let Ω 3 be the family of curves lying in the upper  Fig. 19 Comparison with the exact formula (29) for n = 2 13 half-plane and joining the sets (−1, 0) and (s, ∞). The modulus is invariant under conformal mapping [1,13,24] and hence [7, p. 20] Consequently, the exact value of the constant r is given by where the value of the s + 1 is given by (27). In this paper, the values of the function μ are computed as described in [20].
To test the Algorithm 6.2, we fix θ 1 = 0.5π and θ 3 = 1.5π . Then, we choose values for θ 2 between 0.5001π and 1.4999π . The numerical results obtained with n = 2 13 are shown in Fig. 19. Figure 19 shows the relative error in the computed values of r vs θ 2 (left), the total CPU time (in seconds) required for computing each value of r vs θ 2 (center), and the successive error |r k − r k−1 | for each value of θ 2 vs the iteration number k (right). We see from the figure, with less than 40 iterations, the successive error for the iterative method method is less than 10 −13 for all values of θ 2 except for θ 2 = 0.5001π (red line). As expected, the relative error in the computed values of r is very small when θ 2 is a way from θ 1 and θ 3 . The numerical results obtained with n = 2 10 are shown in Fig. 20.  Fig. 20 Comparison with the exact formula (29) for n = 2 10

The Crowding Phenomenon
According to [6, pp. 20-21], the term crowding was coined in 1980 [15] to describe the error/instability in numerical computing of conformal mapping. Thereafter it has become a "benchmark issue" for all numerical conformal mapping software. As explained in [21, p. 77], mapping a rectangle with aspect ratio m to the unit disk seems to be impossible for m = 24. Problems start already with m = 8 and become more serious with increasing m. The critical value of m depends on the computer floating point arithmetic with 10 < m < 20 for double precision arithmetic [6, pp. 20-21], [21, pp. 75-77].
For rectangle R r in (22), the aspect ratio is m = r for r > 1 and m = 1/r for r < 1. Assume that θ 1 = π/2 and θ 3 = 3π/2 are fixed as above and π/2 < θ 2 < 3π/2. In this subsection, we use the analytic example presented in Sect. 6.4 to find the critical value of r for mapping the quadrilateral (D; 1, i, e iθ 2 , −i) onto the rectangle R r with double precision arithmetic.
In view of (27), we have Thus, by (29) the value of the modulus r can be written in terms of θ 2 , Also, the value of θ 2 can be written in terms of the modulus r , The values of the modulus r obtained with the formula (30) for π/2 + 10 −15 < θ 2 < 3π/2 − 10 −15 are shown in Fig. 21 (left). Similarly, Fig. 21 (right) shows values θ 2 obtained with the formula (31) for 1/12 < r < 12. We see from Fig. 21 (right) that Similarly, for r = 1/12, we have For further investigation, we use the MATLAB symbolic toolbox to obtain more accurate results for the formulas (30) and (31). For r ≥ 1, the values of log 10 (3π/2−θ 2 ) are shown in Fig. 22 (left). It is clear from Fig. 22 (left) that, up to double precision accuracy of the computer, θ 2 = 3π/2 for the value of r as small as r = 13. Similarly, for 0 < r ≤ 1, the values of log 10 (θ 2 − π/2) are shown in Fig. 22 (right). Up to double precision accuracy of the computer, it follows from Fig. 22 (right) that θ 2 = π/2 for r = 1/13. As a consequence, up to double precision accuracy of the computer, mapping a rectangle onto a quadrilateral (D; 1, i, e iθ 2 , −i) is impossible when the aspect ratio of the rectangle is as small as r = 13.
Funding Open Access funding provided by the Qatar National Library.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.