PlgCirMap: A MATLAB toolbox for computing conformal mappings from polygonal multiply connected domains onto circular domains

This paper presents a MATLAB toolbox for computing the conformal mapping from a given polygonal multiply connected domain onto a circular multiply connected domain and its inverse. The toolbox can be used for multiply connected domains with high connectivity and complex geometry. It can be employed also for simply connected domains.


Motivation and significance
Conformal mappings are used to transform two-dimensional domains with complex geometry (physical domains) onto domains with simpler one (canonical domains). Numerous canonical domains have been considered in the literature for conformal mappings of multiply connected domains in the extended complex plane C∪{∞} [15,19,20]. Perhaps the most important canonical domain for simply and multiply connected domains is the circular domain; i.e., a domain all of whose boundaries are circles. This is due to the existence of analytic formulas for several problems in circular multiply connected domains (see the recent monograph [3] and the references cited therein). Furthermore, circular domains are ideal for using Fourier series and FFT [5].
Important examples of complex geometry domains are the polygonal domains, whose boundaries consist of straight line segments. For simply connected domains, the Schwartz-Christoffel (SC, for short) formula provides us with an explicit form of the conformal mapping from the unit disk onto a given polygonal domain [9]. The SC formula for simply connected domains was discovered independently by Christoffel in 1867 and Schwarz in 1869 (see [9, p. 4]). The generalization of this formula to doubly connected domains was due to Komatu [16] in 1945 (see [12, pp. 478-486]). However, the extension of the SC formula to multiply connected domains was established only recently. Indeed, DeLillo, Elcrat and Pfaltzgraff [5] and DeLillo [4] derived SC formulas for conformal mappings from circular domains onto unbounded and bounded polygonal domains, respectively, using the reflection principle. Crowdy [1,2] presented SC formulas for computing such mappings using Schottky-Klein prime functions.
Driscoll [6,7,8] created a MATLAB package called SC Toolbox for computing the conformal mapping from the unit disk onto a given polygonal simply connected domain. The toolbox is a generalization of the Fortran package SCPACK developed by Trefethen [23]. The SC Toolbox has been widely used by many researchers. However, no such toolbox is available so far for polygonal multiply connected domains. The development of such a MATLAB toolbox is the subject of this paper. The proposed toolbox can be used for computing the conformal mapping w = f (z) from a given polygonal multiply connected domain G onto a circular domain D and its inverse z = f −1 (w).

The conformal mapping
Assume G is a given bounded or unbounded polygonal multiply connected domain bordered by m polygons Γ j , j = 1, 2, . . . , m, such that no corner of these polygons is a cusp. For m = 1, the domain G is simply connected. If G is bounded, then we assume that Γ m is the external boundary component and encloses all the other boundary components Γ j , j = 1, 2, . . . , m−1. The total boundary Γ = ∂G = ∪ m j=1 Γ j is oriented such that G is on the left of Γ. Then there exists a conformal mapping w = f (z) from the domain G onto a circular multiply connected domain D bordered by m circles C j , j = 1, 2, . . . , m (see [12, p. pp. 488-496] and [25, pp. 118-127]). The mapping f extends to the boundary of G with f (Γ j ) = C j , j = 1, 2, . . . , m. We assume that the circular domain D is bounded when G is bounded and D is unbounded when G is unbounded. The total boundary C = ∂D = ∪ m j=1 C j has the same orientation as Γ.
When G is bounded, the conformal mapping f is uniquely determined by assuming that the external boundary C m = f (Γ m ) of D is the unit circle and where α is a given point in the domain G. The condition (1) can be replaced with the condition where β is a given point on the external boundary Γ m . On the other hand, if G is unbounded, then the mapping function f is uniquely determined by assuming that near infinity. Alternatively, f is uniquely determined by assuming that C m = f (Γ m ) is the unit circle and where β is a given point on the boundary Γ m . The SC formulas derived in [1,2,4,5] can be used to compute the inverse mapping z = f −1 (w) from the circular domain D onto the polygonal domain G.
However, using these SC formulas requires solving a system of non-linear equations to determine the preimages of the polygons' vertices as well as the centers and the radius of the circles. Solving such a nonlinear system of equations is still a challenging problem. On the other hand, conformal mappings from multiply connected domains onto circular domains can be computed using Koebe's iterative method [14] (see [12, § 17.7]). As a special case, Koebe's method can be used to compute the conformal mapping w = f (z) from a given polygonal multiply connected domain G onto a circular multiply connected domain D. A fast implementation of Koebe's iterative method using the boundary integral equation with the generalized Neumann kernel is given in [22] (see also [21]). The method presented in [22] can be used also to compute the inverse map z = f −1 (w) from the circular domain D onto the polygonal domain G. Further, this method can be applied even for domains with high connectivity and complex geometry, see [21,22]. More recently, another efficient numerical method for computing the conformal mapping from circular domains onto polygonal domains based on rational approximations has been presented by Gopal and Trefethen [10] and Trefethen [24].
In the presented PlgCirMap toolbox, the above described conformal mapping w = f (z) from the polygonal domain G onto the circular domain D as well as its inverse z = f −1 (w) from D onto G will be computed using the boundary integral method presented in [22].

Software Architecture
The PlgCirMap is a MATLAB toolbox that consists of different MATLAB functions. The main functions in this toolbox are plgcirmap, evalu, evalud, and plotmap. The function plgcirmap itself depends on three main functions, mainmap, cirmapb, and cirmapu (see Figure 1). The inputs for the function plgcirmap are a cell array ver containing the vertices of the polygons and a point alpha in the polygonal domain G. The default values of the parameters for the numerical calculations are set in the function plgcirmap. The numerical implementation of Koebe's iterative method is given in the function mainmap. Koebe's iterative method requires computing conformal mappings from bounded simply connected domains onto the unit disk and computing conformal mappings from unbounded simply connected domains onto the exterior of the unit disk (see [22] for more details). Such conformal mappings will be computed using the functions cirmapb and cirmapu, respectively. The output of the plgcirmap will be a MATLAB object f containing the required information about the conformal mapping f and its inverse f −1 . From the object f, we can compute the values of the conformal mapping and its inverse using the function evalu. The values of the first derivatives of f and f −1 can be computed using the function evalud. Finally, the function plotmap uses the object f to visualize the conformal mapping f and its inverse f −1 .

Software Functionalities
The PlgCirMap toolbox is used to compute and visualize the conformal mapping w = f (z) from a given polygonal multiply connected domain G onto a circular multiply connected domain D and its inverse z = f −1 (w). All conditions (1) can be implemented in the toolbox. PlgCirMap can be used also for polygonal simply connected domains. To use the PlgCirMap toolbox, the boundary of the polygonal domain is assumed to have no cusps or slits.

Parameters' default values
In the PlgCirMap toolbox, the mapping function w = f (z) and its inverse z = f −1 (w) are computed using the boundary integral method presented in [22]. The method is based on a fast numerical implementation of Koebe's iterative method using the boundary integral equation with the generalized Neumann kernel [21,22]. Assume that each polygon Γ j has ℓ j ≥ 3 vertices. We discretize each segment of the polygon Γ j with n graded mesh points so that Γ j is discretized by ℓ j n graded mesh points (see [17,18] for details on how the graded mesh points are chosen). Then, applying the Nyström method with the trapezoidal rule reduces the integral equation to a linear system which is solved iteratively by the MATLAB function gmres. The matrix-vector product in the GMRES method is computed using the MATLAB function zfmm2dpart from the MATLAB toolbox FMMLIB2D [11]. The computational cost of the method is O(m 2 ℓn + mℓn log n) where ℓ = max 1≤j≤m ℓ j (see [21,22] for details).
In the PlgCirMap toolbox, the default values of the parameters for the numerical calculations are set in the function plgcirmap as follows.
1. The default value of n, the number of discretization points in each side of the polygons, is set to n=2^9. In fact, accurate results can be obtained even for the values of n as small value as n=2^5. Increasing the value of n leads to increased accuracy of the obtained results. However, choosing very large values of n should be avoided since it could cause a problem with the convergence of the FMM functions in the toolbox FMMLIB2D [11].
2. For the FMM function zfmm2dpart, we set the default value iprec=4, which means the accuracy of the FMM is 0.5 × 10 −12 . The accuracy of the obtained results can be improved by choosing iprec=5 (the accuracy of the FMM will be 0.5 × 10 −15 ). However, there may be a problem with the convergence of the FMM if we choose iprec=5, especially when n is too large.
3. In the MATLAB function gmres, the default tolerances of the GMRES method is set to gmrestol=0.5e-12 and the default maximum number of iterations allowed is set to gmresmaxit=100. The GMRES is used without restart.
4. For Koebe's iterative method, the default tolerance and the default maximum number of iterations allowed are set to koebetol=1e-12 and koebemaxit=100, respectively. Remarks.
1. The default value for maximum number of iterations allowed for both the GMRES method and Koebe's iterative method is 100. However, for several numerical experiments with the PlgCirMap toolbox, both methods converges with less than 100 iterations.
2. If we choose very large values of n and/or iprec=5, the FMM functions in the toolbox FMMLIB2D might cause the computer to crash. Unfortunately, no warning message will be displayed and sometimes one need to restart MAT-LAB. In such case, we need to reduce the value of n. To compute the conformal mapping with the normalizations (2) or (4), we call the function plgcirmap as: which means that the vertex k on the polygon Γ m will be mapped to 1. The object f, which is a MATLAB struct with several fields, contains the data of the conformal mapping w = f (z) from G onto D as well as its inverse z = f −1 (w) from D onto G.

The function evalud
The function evalud is used to compute the values of the first derivatives of the mapping function f and its inverse f −1 . The values of f ′ (z) at a vector of points z in G can be computed via evalud(f,z,'d') and the values of (f −1 ) ′ (w) at a vector of points w in D can be computed by evalud(f,w,'v').

The function plotmap
The function plotmap is used to visualize the conformal mapping f and its inverse f −1 . To plot rectangular grids in the polygonal domain G and their images in the circular domain D, we call plotmap(f,'d','rec',n1,n2) where n1 is the number of horizontal lines and n2 is the number of vertical lines. We can also plot rectangular grids in D and their images in G through plotmap(f,'v','rec',n1,n2). Similarly, polar grids can be plotted via plotmap(f,'d','plr',n1,n2) and plotmap(f,'v','plr',n1,n2) where n1 is the number of circles and n2 is the number of rays. To plot the domains G and D without grid-points, we call plotmap(f).

Comparison with Schwarz-Christoffel Toolbox [8]
The PlgCirMap toolbox can be used for computing the conformal mapping from a given polygonal simply connected domain G onto the unit disk (for bounded G) or the exterior of the unit disk (for unbounded G). For such case, the accuracy of the PlgCirMap toolbox can be compared against the well known SC Toolbox [8].
Consider the simply connected domain G interior to the polygon with the vertices 1.5i, −1 + 1.5i, −1 − i, 1.5 − i, 1.5, and 1 (see Figure 2). The SC Toolbox can be used to map the unit disk onto the domain G such that 1 on the unit circle is mapped to the last vertex of the polygonal which is also 1. So, we shall assume here that f is the conformal mapping from the polygonal domain G onto the unit disk D normalized by (2)   The resulting figure is shown in Figure 2 (center).
Next, we use the SC Toolbox to compute the inverse mapping f −1 from the unit disk D onto the polygonal domain G. First, we set the accuracy of the SC toolbox to 10 −14 by calling >> options = scmapopt('Tolerance',1e-14); Then, we use the SC toolbox to compute a MATLAB object fisc by calling >> p=polygon(ver{1}); >> fisc = diskmap(p,options); >> fisc = center(fisc,alpha); The plot of the image of the orthogonal polar grids in the disk D under the inverse mapping f −1 computed by the SC toolbox is shown in Figure 2  Then, we compute the maximum norm between the computed values as >> E_1 = norm(prevertsc-prevertpcm,inf) The obtained maximum norm is E 1 = 1.0579 × 10 −12 .
Next, we choose a set of points >> zz = 0.6.*exp(i.*linspace(0,2*pi,1000)); in the polygonal domain G. To compute the values of the conformal mapping f at these points zz using the SC toolbox, we call the function evalinv(fisc,zz). For the PlgCirMap toolbox, the values of the mapping f at the points zz are computed by calling evalu(f,zz,'d'). Then, we compute the maximum norm between the computed values as >> E_2 = norm(evalinv(fisc,zz)-evalu(f,zz,'d'),inf) The obtained maximum norm is E 2 = 5.3952 × 10 −13 .
We also compare the two toolboxes by choosing a set of points >> ww = 0.9.*exp(i.*linspace(0,2*pi,1000)); in the unit disk D. Then, we compute the maximum norm between the values of the inverse map f −1 computed at the points ww by the two toolboxes as The obtained maximum norm is E 3 = 1.0291 × 10 −12 . Finally, we check the accuracy of each toolbox separately. We use both toolboxes to compute approximate values for f −1 (f (zz)) and f (f −1 (ww)). Then, we compute the maximum error norm in the computed values for the SC toolbox by As we see from this example, there is a very good agreement between the results obtained by the well developed SC toolbox and the presented PlgCirMap toolbox (for the default value of n, n = 2 9 ). Further, we will have a good agreement between the results obtained by the two toolboxes even for small values of n. Indeed, if we change the value of n in the MATLAB function plgcirmap to n = 2 5 , then the above computed norms will be as follows: E 1 = 9.6102 × 10 −8 , E 2 = 1.0283 × 10 −8 , E 3 = 6.3408 × 10 −8 , E 4 = 4.7186 × 10 −9 , and E 5 = 8.2561 × 10 −9 .

Illustrative Examples
In this section, we use the PlgCirMap toolbox to compute the conformal mappings for two domains. More examples for both bounded and unbounded domains are available in https://github.com/mmsnasser/PlgCirMap.

A bounded multiply connected domain
For the first example, we consider a bounded polygonal multiply connected domain of connectivity 17. In the MATLAB code given below, we first define the vertices of the polygons and we choose a point α in the domain G. Then, we call the function plgcirmap to compute the conformal mapping with the normalization (1). The function plotmap is then called to visualize the conformal mapping f and its inverse f −1 as in Figure 3.  20,20); plotmap(f,'d','plr', 15,25); plotmap(f,'v','plr',10,20);

An unbounded multiply connected domain
In the second example, we consider an unbounded polygonal multiply connected domain of connectivity 24. The MATLAB code for this example is given below where the normalization (3) is used. Figure 4 shows the obtained figures.

Impact
Conformal mappings are a powerful tool to solve several problems in the fields of science and engineering involving the Laplace equation due to its invariant under conformal mappings. With the help of conformal mappings, solving the Laplace equation in domains with complex geometry (physical domains) is reduced to solving this equation in domains with simpler geometry (canonical domains). The simple geometry of the circular domain makes it an important canonical domain from both physical and computational points of view. For example, D. Crowdy with several collaborators have recently presented analytic formulas for several problems of fluid mechanics in circular multiply connected domains (see e.g., [3] and the references cited therein). With the help of the presented toolbox, such analytic formulas can be used also for polygonal domains (we refer to [13] for an example of such applications of the toolbox).
7 Final remarks and suggestions for future improvements 1. The numerical method used in the PlgCirMap toolbox to compute the conformal mapping is based on the boundary integral method presented in [22]. Thus, the accuracy of the toolbox depend on the accuracy of the numerical solution of the used boundary integral equation. In the current version of the toolbox, the integral equation is solved using the Nyström method with the trapezoidal rule based on graded mesh points (see [17] for details). Improving the accuracy of the numerical solution of the integral equation will improve the accuracy of the presented toolbox. This could be considered in the future.
2. The method presented in [22] can be used for general multiply connected domains with smooth or piecewise smooth boundaries. As a result, the PlgCirMap toolbox can be generalized to multiply connected domains other than polygonal domains. In particular, this toolbox can be generalized easily to compute the conformal mapping from multiply connected domains with circular-arc boundaries onto circular multiply connected domains.
3. The accuracy of the numerical methods for computing the conformal mapping from a circular domain D onto an elongated domain G is seriously affected by what is known as the crowding phenomenon (see [9, §2.6]). Although the method used in this toolbox is based on computing the conformal mapping from the elongated domain G onto the circular domain D, it still affected by crowding. One possible way to improve the accuracy of the presented toolbox for elongated domains in future is to use the domain decomposition method or to consider canonical domains other than the circular domain (see [9]).