Log In Sign Up

Numerical conformal mapping with rational functions

by   Lloyd N. Trefethen, et al.

New algorithms are presented for numerical conformal mapping based on rational approximations and the solution of Dirichlet problems by least-squares fitting on the boundary. The methods are targeted at regions with corners, where the Dirichlet problem is solved by the "lightning Laplace solver" with poles exponentially clustered near each singularity. For polygons and circular polygons, further simplifications are possible.


page 1

page 2

page 3

page 4


Solving Laplace problems with corner singularities via rational functions

A new method is introduced for solving Laplace problems on 2D regions wi...

AAA-least squares rational approximation and solution of Laplace problems

A two-step method for solving planar Laplace problems via rational appro...

Lightning Stokes solver

Gopal and Trefethen recently introduced "lightning solvers" for the 2D L...

Exponential node clustering at singularities for rational approximation, quadrature, and PDEs

Rational approximations of functions with singularities can converge at ...

Model order reduction of layered waveguides via rational Krylov fitting

Rational approximation recently emerged as an efficient numerical tool f...

Solving Laplace problems with the AAA algorithm

We present a novel application of the recently developed AAA algorithm t...

A linear barycentric rational interpolant on starlike domains

When an approximant is accurate on the interval, it is only natural to t...

1 Introduction

The association of conformal mapping with Laplace problems goes back to the beginning of the subject henrici ; schiffer . When this connection has been exploited for the numerical computation of conformal maps, it has usually been by way of integral equations such as those of Berrut, Gerschgorin, Lichtenstein, Symm, and Warschawski gaier ; henrici ; tref ; wegmann . However, almost all domains of interest in practical conformal mapping contain corners, and for such domains, it has recently been discovered that rational functions with exponentially clustered poles offer a powerful alternative approach to solving Laplace problems lightning . The aim of this paper is to bring these ideas together to present simple new algorithms for numerical conformal mapping. A second, different aspect of rational functions introduced in gt also comes into play: the compression of a computed conformal map to a more compact form by means of the AAA algorithm AAA ; lawson .

As shown in lightning , rational functions can achieve root-exponential convergence for almost arbitrary corner singularities, i.e., as a function of the rational degree  for a constant , and no analysis of the singularities is needed. This means that our method does not care if the boundary components are straight or curved, so long as they are smooth. Unbounded domains are also not a problem.

For simply-connected domains with smooth boundaries, the idea of conformal mapping via Laplace problems can be realized by means of polynomials instead of rational functions, and we begin in Section 2 by presenting an algorithm of this kind. Though this method is sometimes effective, it often fails for reasons related to the “crowding phenomenon” of numerical conformal mapping. Methods based on integral equations are superior, and such a method is the default used in the new Chebfun code conformal. The picture changes sharply when we turn to rational functions in Section 3. We show that a conformal map of a simply-connected domain with corners can be computed in a few lines of code by calling the MATLAB “lightning Laplace solver” laplace lightcode .

Both the polynomial and rational function methods are readily generalized to doubly connected domains, and these generalizations are presented for smooth boundaries in Section 4 and for boundaries with corners in Section 5. Domains of higher connectivity are considered in Section 6, and numerical details of our methods are discussed in Section 7. Section 8 discusses the further simplifications available for polygons and circular polygons if one uses a modified method based on explicit information about corner singularities. Conclusions are presented in Section 9.

2 Smooth simply connected domains

Let be a simply connected domain bounded by a Jordan curve , and let and be the open unit disk and the unit circle. Without loss of generality we assume and seek the unique conformal map such that and . By the Osgood–Carathéodory theorem henrici , extends to a homeomorphism of onto . It follows that is a nonzero analytic function on that is continuous on and has real part for and imaginary part at . If we write as , where and are real harmonic functions in , continuous on , then is the unique solution of the Dirichlet problem


and is its unique harmonic conjugate in with . Combining these elements, we see that takes the form


Note that is the Green’s function of with respect to the point , so is nothing more than the exponential of the analytic extension of the Green’s function. See (henrici, , Thm. 16.5a) and (schiffer, , p. 253).

To compute the conformal map , therefore, it is enough to solve the Dirichlet problem (1) for and then find its harmonic conjugate . In this paper we approximate solutions to Dirichlet problems by the real parts of analytic functions defined by series involving polynomials, rational functions, or other elementary terms. The harmonic conjugates come for free from the same expansion coefficients.

In this section we consider the simplest case, in which is well-behaved enough that it is effective to work with polynomials. We first outline the algorithm in the monomial basis and then mention the variation involving Arnoldi factorization that our codes actually employ. Following series , let be approximated as


for some degree and complex coefficients , which will depend on . Writing , we have


With , this gives us real unknowns and . These are determined numerically by sampling the function in points on and solving an matrix least-squares problem of the form ; we take . The harmonic conjugate is then approximated by the series


and comes from (2). Note that if is analytic, then can be analytically continued to a neighborhood of , and geometric convergence of (4) and (5) is guaranteed (in exact arithmetic) as  walsh , assuming the boundary sampling is fine enough that the discrete least-squares problem on approximates the continuous one. (There are also special cases in which is not smooth and yet is still analytic everywhere on , such as a rectangle, an equilateral triangle, or the “circular L-shape” consisting of a square from which a quarter-circular bite has been removed.)

By the method just described, if is well enough behaved, one can compute a polynomial-based approximation of the conformal map . The degree will usually have to be high. However, as described in gt , the computed map can be compressed to a much more compact form as a rational function by means of the AAA approximation algorithm AAA . Moreover, the AAA method can be invoked a second time to compute a rational approximation of the inverse map . As emphasized in gt , when AAA approximation is brought into play, the traditional fundamental difference between conformal maps in the forward and inverse directions ceases to be very important. Once one is available, the other is readily computed too.

The presentation of the last page has been framed with respect to monomial bases Unless is close to a circle about the origin, however, the monomials will be very ill-conditioned, and it is better to use an approximation to a family of polynomials orthogonal over . This we do by carrying out an Arnoldi iteration as described in VmeetsA . Mathematically the two approaches are equivalent.

We summarize the algorithm, which we call Algorithm P1 (P for polynomial, 1 for simply connected), as follows. This is a high-level description, omitting many details, and the parameters indicated are engineering choices, not claimed to be optimal. For reasons discussed in gt , high accuracy tends to be problematic in numerical conformal mapping, and we work by default to a tolerance of .

Algorithm P1.

For rounded values of until convergence or failure,

- Set and sample in points of ;

- Solve an least-squares problem for the approximation ;

- For stability, do this via Arnoldi factorization as described in VmeetsA .

Form from and from ;

Approximate by a AAA rational approximation;

Approximate by another AAA rational approximation.

This algorithm is available in Chebfun chebfun in the command conformal, introduced in October, 2019, if it is called with the flag ’poly’. It is not the default algorithm used by conformal

, as we shall explain in a moment.

Figure 1 illustrates Algorithm P1 applied to a smooth domain with five-fold symmetry, showing the output from the Chebfun commands

     C = chebfun(’exp(1i*pi*t)*(1+.15*cos(5*pi*t))’, ’trig’);
     [f, finv] = conformal(C, ’poly’, ’plots’);

The objects f and finv are function handles corresponding to rational approximations of degrees 26 and 22, respectively, of the functions and . These degrees are remarkably low. The computation takes about seconds on a desktop machine, plus another second for plotting, and the maximal error on the boundary is about .

Figure 1: Conformal map of a smooth domain computed to about -digit accuracy in seconds on a desktop machine by Algorithm P1 via the Chebfun conformal(...,’poly’) command. The red dots mark poles of the rational representations of and : just poles for and for . The closest poles lie at distance from on the left and at distance from on the right.
Figure 2: If the lobes are made deeper by changing the coefficient of the cosine term from to , Algorithm P1 fails because it would require a polynomial of degree in the tens of thousands. Instead this figure was computed via the Kerzman–Stein integral equation, the default algorithm of Chebfun conformal. The closest poles now lie at distances from on the left and from on the right. The AAA rational representations are still remarkably economical.

All this sounds good, but in fact, Algorithm P1 is limited in its capabilities, and will fail for many smooth domains. For example, it fails if the indentations in the domain of Figure 1 are made twice as deep by changing the coefficient to . The problem in such a case is that the polynomial degree would have to be in the thousands, a consequence of the large distortions in conformal mapping that go by the name of the crowding phenomenon wegmann . If is increased further to , we get the domain shown in Figure 2, and this was certainly not mapped with Algorithm P1. The default method of conformal was applied instead, involving a discretization of the Kerzman–Stein integral equation ks ; kt followed by AAA rational compression in both directions. For this domain, Algorithm P1 would require a polynomial of degree in the tens of thousands.

There is a literature on the crowding phenomenon, and in particular, see Theorems 2–5 of gt . These establish that if contains an outward-pointing finger of length-to-width ratio as defined in a certain precise sense, then must be at least as large as order at some points of , with singularities (or at least points of non-univalency) at distances no greater than order from , and that any useful polynomial approximation must have degree at least on the order of . (In Figure 2 we have and .) Analogous if somewhat weaker effects occur for and its polynomial approximations if there are inward-pointing fingers. The geometric convergence mentioned earlier still applies, but is useless in practice because the constant exceeds by only a small amount. These effects are reflected in the positions of the red dots marking poles in Figures 1 and 2, highlighting the profound difference between polynomial and rational approximations. In the right image of each pair, the poles cluster near points of that map to outward-pointing lobes of . Similarly, the poles in the left image cluster near points of at inward-pointing fingers.

We emphasize that the effects just discussed will cause trouble for any method based on polynomial approximations, no matter how stably implemented, for they have nothing to do with the choice of polynomial basis. Our Arnoldi factorization makes the implementation quite stable, but still it will fail for domains whose maps involve large distortions.

Variants of Algorithm P1 can be used to solve related problems involving exterior domains. To map the interior of to the exterior of , we change the factor to in (2); the Dirichlet problem (1) remains unchanged. To map the exterior of to the exterior of , we replace positive by negative powers in (3)–(5). To map the exterior of to the interior of , we make both of these changes.

We now turn now to the much more powerful Algorithm R1 based on rational functions.

3 Simply connected domains with corners

Let the boundary of  be a polygon or circular polygon or other Jordan curve consisting of a finite set of analytic curves meeting at corners . (A circular polygon is a finite chain of circular arcs.) As before, we assume and seek the conformal map with and .

Equations (1) and (2) still hold: we need to solve a Dirichlet problem, but now, there will be singularities at the corners. Instead of (3), we consider


with complex coefficients and ; again, the polynomial part of the sum will be implemented in practice by an Arnoldi factorization VmeetsA . Since the publication of a paper by Newman in 1964 newman , it has been recognized that rational functions of the form of the first sum in brackets can approximate certain singularities with root-exponential convergence, if the poles  are clustered exponentially near the singular points with a spacing that scales with . In lightning it was proved that this effect applies generally to solutions of Lapalce Dirichlet problems defined by piecewise analytic functions meeting at corners. In particular, Theorem 2.3 of lightning implies that approximations of the form (6) can converge root-exponentially for the conformal map . (The proof given in lightning assumes is convex, but it is believed that the result holds generally.) Since the appearance of lightning , a MATLAB “lightning Laplace solver” code laplace has been developed to solve Dirichlet problems by this method lightcode .

We summarize our algorithm, Algorithm R1.

Algorithm R1.

Solve problem  by the lightning method and form from (2);

Approximate by a AAA rational approximation;

Approximate by another AAA rational approximation.

Using laplace, all this can be done in four lines of MATLAB, assuming a specification of and a tolerance tol have been given. If is just a polygon, then P

can be just a vector of vertices. By default we take


     [~, ~, w, Z] = laplace(P, @(z) -log(abs(z)), ’tol’, tol);
     W = Z.*exp(w(Z));
     f = aaa(W, Z, ’tol’, tol);
     finv = aaa(Z, W, ’tol’, tol);

This is the central message of this paper: the conformal map of a domain with corners can be computed in this simple fashion. Section 7 will discuss certain numerical details.

Figure 3: Conformal map of a pentagon computed in about secs. by Algorithm R1.

Figure 3 illustrates Algorithm R1 applied to a regular pentagon, a polygonal analogue of the domain of Figure 1. Five-digit accuracy is achieved in about secs. of desktop time, with poles in both directions clustering near the five singular points. There are four poles near the vertex at , for example, whose distances from the vertex are approximately

     0.0097,  0.044,  0.14,  0.49.

Though these poles were chosen adaptively by the AAA algorithm, based on no mathematical analysis, they show exponential clustering with a ratio of successive distances on the order of .

Figure 4: Mapping of an L-shaped domain by Algorithm R1. The map is nonsingular at the salient corners, so on the left, poles cluster only at the reentrant corner, whereas on the right they cluster at six points.

Figure 5: The rational approximations produced by Algorithm R1 can be evaluated in just microseconds. This figure shows images under

of 10,000 uniformly distributed random points in the unit disk, computed in about

secs. total. The nonuniform shading illustrates the exponential distortion effects in conformal maps, discussed in Section .

Figure 4 shows results of an analogous computation, taking about half a second, for an L-shaped domain. On the left, poles cluster near the reentrant corner at distances approximately

    0.000053, 0.00030, 0.00069, 0.0071, 0.0013, 0.0022, 0.0035,
    0.0054, 0.081, 0.012, 0.017, 0.024, 0.034, 0.047, 0.066,
    0.091, 0.12, 0.17, 0.23, 0.31, 0.42, 0.56, 0.76, 1.04, 1.66.

These show exponential spacing with a factor of about 1.4. The reduction of the factor from 3 to 1.4 as the number of poles near a corner increases from 4 to 25 is consistent with root-exponential clustering theorems of lightning .

We emphasize the extraordinary economy of the representations of conformal maps computed by Algorithm R1, since they are nothing more than rational functions of a degree typically less than 100.  Figure 5 shows images of points in the unit disk, computed in secs. all together.

Figure 6: Conformal map of a circular polygon computed by Algorithm R1. This shape is the first example presented by Howell in howell .

Figure 7: Another map of a circular polygon computed by Algorithm R1. Each boundary component is a circular arc of radius .

These examples involve polygons, but Algorithm R1 is equally applicable to other domains bounded by analytic arcs meeting at corners. For example, a circular polygon is a generalization of a polygon in which each side may be a straight segment or a circular arc. Whereas polygons are conventionally mapped by the Schwarz–Christoffel integral, circular polygons can be mapped by the so-called Schwarzian differential equation SCbook . This idea goes goes back to Schwarz himself 150 years ago, but has rarely been attempted numerically; the two realizations we know of are bg and howell . Figures 6 and 7 show maps of circular polygons computed by the completely different approach of Algorithm R1, each in less than a second on a desktop. The laplace code used for these computations includes a convenient syntax for specifying circular boundary arcs via their initial point and radius of curvature. Thus the mapping of Figure 7 was computed by the four-line command sequence listed above preceded by the specification

     P = { [1 -2], [1i -2], [-1 -2], [-1i -2] };

Figure 8: Exterior map for the region of Figure 7 computed by a variant of Algorithm R1.  The poles of (6) now lie interior to the boundary curve.

Just as noted at the end of the last section for smooth domains, variants of Algorithm R1 can be developed for exterior geometries. To map the interior of to the exterior of , we change the factor to in (2) without altering the Dirichlet problem (1). To map the exterior of to the exterior of , we place poles in (6) interior to rather than exterior and also change to . To map the exterior of to the interior of , we make both of these changes. Figure 8 illustrates an exterior-exterior map for the same curve as in Figure 7. (Composing interior and exterior maps for a given curve leads to the idea of “conformal welding” bishop ).)

4 Smooth doubly connected domains

Conformal mapping via Dirichlet problems, like the idea of Green’s functions it is based upon, does not just apply for simply connected domains. Following the pattern of the last few pages, we will first introduce the mapping of smooth doubly-connected domains in two alternative geometries, then generalize in the next section to doubly-connected domains with corners.

Suppose is the annulus bounded by an outer Jordan curve and an inner Jordan curve , both enclosing the origin . Let be a conformal map of onto the circular annulus with and for some , known as the conformal modulus

. Such a map is unique up to one real degree of freedom corresponding to a rotation. The value of

is unknown a priori and will be determined as part of the solution. Writing again in the form (2) shows that we must find a harmonic function in satisfying

see (schiffer, , p. 253ff).

Now for any value of , not just the one we are looking for, there will be a unique solution to the Dirichlet problem (4axler . If and are smooth, we could approximate this solution by functions of the form


for some real , with convergence as (see the “Logarithmic Conjugation Theorem” of axler , and the final section of the same paper for this application to doubly connected conformal mapping). However, if is not an integer, then the analytic function of (2) associated with this formula is not single-valued. Even if is an integer, will fail to be one-to-one in if . Only with do we get the conformal map we seek. And thus we find that the problem to be solved is to find functions of the form


converging to the boundary values (4) for some as , and the value of for which this is possible is the conformal modulus of . The analogues of (4) and (5) will be


with .

Figure 9: Above, map of a smooth annulus computed by Algorithm P2. Poles of the rational representations of and appear both inside and outside the annuli. Below, an unbounded smooth doubly-connected geometry mapped via the variant of Algorithm P2. The point is now in the interior of , and its image on the right is a simple pole of .

To make (4) and (9) the basis of an algorithm, we set up a least-squares problem whose unknowns are the coefficients and and also the conformal modulus . The number of and coefficients is , and the total number of unknowns is . The least-squares matrix is of dimensions , with the rows corresponding to the sample points on the boundary curves and the first columns containing the values of and at these points. The final column consists of zeros and ones, corresponding to the appearance of at one boundary curve but not the other in (4), and the right-hand side is the -vector consisting of samples on the boundary of . Here is the algorithm:

Algorithm P2.

For rounded values of until convergence or failure,

- Set and sample in points each of and ;

- Solve a least-squares problem for the coefficients and of and the number ;

Form and from and from ;

Approximate by a AAA rational approximation;

Approximate by another AAA rational approximation.

Chebfun’s conformal command executes this algorithm when it is called with two boundary arcs rather than one,111[Note to referees: this extension is not yet in the publicly available version of Chebfun but will be introduced shortly.] and the upper half of Figure 9 shows the result computed in about sec. by these Chebfun commands:

    z = chebfun(’exp(1i*pi*z)’,’trig’);
    C1 = z*abs(1+.1*z^4); C2 = .5*z*abs(1+.2*z^3);
    conformal(C1, C2, ’plots’);

So far as I know, Figure 9 is the first illustration ever published of rational approximations with free poles used to represent conformal maps of doubly-connected domains.

There is an alternative doubly-connected geometry which, while mathematically equivalent to the one we have just treated, may in practice be worth considering separately. Suppose that , rather than being the annular region between two concentric curves and , is the infinite region of the extended complex plane exterior to two disjoint curves and . If encloses and encloses , then we can write the conformal map of onto a circular annulus in the form


where is harmonic in and satisfies, in analogy to (4),

for a value of that again will be uniquely determined. To fashion this into an elementary algorithm based on series, we replace the Laurent polynomials of (9) by polynomials involving negative powers of and . The algorithm is otherwise the same, and the second row of Figure 9 shows a computed example. Note that , , and are regular at , where they take the values , , and , respectively. Thus the inverse map has a pole in the circular annulus at the point , as appears in the plot.

As with its simply connected analogue P1 of Section 2, Algorithm P2 is far from robust, and it will quickly fail if applied to regions much more complicated than those of Figure 9. Superior performance can undoubtedly be achieved with integral equations, though this has not been implemented in Chebfun.

5 Doubly connected domains with corners

By a doubly connected domain with corners, we mean the region lying between two Jordan curves and as in the last section, where now, each curve is a finite chain of analytic arcs. To conformally map such a domain, we solve the Dirichlet problem of section 4 by the lightning Laplace solver of section 3, placing poles both outside the corners of and inside the corners of . The algorithm can be summarized as follows.

Algorithm R2.

Solve by the lightning method and form from (2);

Approximate by a AAA rational approximation;

Approximate by another AAA rational approximation.

The upper part of Figure 10 illustrates such a map for a region that is a polygonal analogue of the region of the first map of Figure 9. As in the last section, there is a straightforward analogue of Algorithm R2 for an unbounded region exterior to two disjoint curves with corners, and this is illustrated in the lower part of Figure 10. Such domains may be of particular interest for the calculation of constants associated with low-rank separability, as discussed by Beckermann, Townsend, and Wilber bt ; tw .

Figure 10: Doubly connected domains with corners mapped by Algorithm R2 and its variant. The function in the lower pair again has a simple pole in the domain, shown by a dot.

For polygonal boundaries as in Figure 10, there is a doubly-connected Schwarz–Christoffel formula dep ; SCbook ; henrici which has been implemented numerically by Däppen daeppen and Hu hu , the latter with a Fortran package DSCPACK.

6 Multiply connected domains

Domains of connectivity higher than two can also be conformally mapped. The Green’s function idea of (4) generalizes to a formula with different constants on each boundary component, as described in henrici and schiffer . Various numerical methods have been developed, and some of them are reviewed in bds . The Schwarz–Christoffel formula has generalizations that have been developed by DeLillo, Crowdy, and others: crowdy and dep2 are just two of many papers in this area. Recently Nasser has released a code PlgCirMap for computing maps of multiply connected domains by means of a fast implementation of an iterative method due to Koebe nasser .

A difficulty with multiply connected conformal mapping is the lack of an entirely natural target domain. In the doubly connected case, the reduction to an annulus gives a useful interpretation in almost any application. For higher connectivity, standard choices of target domains are regions with holes in the form of circular slits (as delivered by the Green’s function), radial slits, or disks. Physical interpretations in such regions are more strained, however, and it is not always clear that much is gained by transplanting a problem such as a partial differential equation from its native multiply-connected geometry to one of these alternatives.

Figure 11: A multiply connected example domain of Nasser nasser . The map is first computed with his PlgCirMap code, and the forward and inverse maps are then approximated by rational functions, with poles as shown. The rational approximations speed up evaluations by a factor of .

In view of these considerations, we have not extended our numerical methods based on rational functions to multiply-connected regions. If one has computed such a map by an existing method, however, then the possibility remains of compressing it by a subsequent rational approximation as proposed in gt . Figure 11 shows an example. A 4-connected conformal map has been computed by PlgCirMap, and AAA approximation compresses the result to a much more efficient form. Typical times for evaluation of the mapping function or its inverse at each point for this example are on the order of 200 microseconds for PlgCirMap and 10 microseconds for the rational approximations.

7 Numerical details

We now mention some numerical details associated with the algorithms presented in the previous sections. Our implementations are for the most part experimental, however, not claimed to be close to optimal, and we trust that improvements in robustness and accuracy will come in the future.

Accurate data and fine sampling on the boundary. When failure occurs in a lightning Laplace solution, whether for conformal mapping or more generally, the reason is usually that the sampling on the boundary is not fine enough or the sampled data are not accurate enough. In particular, when poles are exponentially clustered near a vertex, it is crucial that the sample points be clustered on the same scale or finer. The laplace software does this automatically according to a scheme described in lightning , and we used the same strategy for the doubly-connected examples of section 5. If the sample points are not dense enough, two things may go wrong: the Laplace problem may be solved inaccurately; and, even if it is accurate, the AAA approximation may fail, with spurious poles appearing in regions where one expects analyticity gt

. (By spurious poles we mean poles of very small residue, also known as “Froissart doublets” since they are paired with nearby zeros that nearly cancel their effect.) The latter problem also arises when the grid is fine but the data are insufficiently accurate relative to the scale of the approximation error, so that the AAA algorithm is effectively attempting to approximate a nonanalytic function. Conversely, if the sample points are dense enough and the data accurate enough, spurious poles rarely appear. The mechanism for this is that AAA algorithm delivers a rational approximation with a certain kind of (linearized) optimality. A spurious pole cannot improve an approximation much, so it is unlikely to be present in an optimal approximation. Success cannot be guaranteed, however, as is well known in the theoretical literature on spurious poles in (mostly Padé) rational approximations 

stahl .

Parameters of AAA approximation. As described in AAA and lawson , Chebfun’s aaa code for AAA approximation has some adjustable parameters. The overall tolerance can be specified, and we have set it equal to the tolerance fixed for the lightning Laplace solver. A “cleanup” option to remove spurious poles is invoked by default, and we have left this in play with its tolerance set equal to the overall AAA tolerance. There is also an option to improve the approximation closer to minimax by a Lawson iteration, and this we have not invoked, because experiments suggest it is ineffective for these problems. All of these choices are experimental, and we have not attempted an analysis that might justify or modify them. As mentioned in the last paragraph, AAA approximation is far from bulletproof. We hope that in the next few years improvements to AAA will be developed, which would have immediate impact on the robustness of our conformal mapping algorithms.

Accuracy and rate of convergence: theory and practice. The theorems of lightning guarantee that Algorithm R1 must converge root-exponentially, i.e. at a rate for some , so long as is bounded by a finite collection of analytic curves meeting at corners. (We note as before that the statements in lightning assume is convex, but it is believed that this condition is not actually necessary.) Sometimes, however, convergence stagnates at a low accuracy such as six digits. One reason is that singularities at corners cause ill-conditioning of the conformal map, as does the crowding phenomenon. (A corner of interior angle with will limit the achievable number of digits to on the order of times that of machine precision.) Another is that our method relies on the solution of least-squares problems defined by ill-conditioned matrices, though this does not seem to cause much trouble in practice, a phenomenon analyzed in the context of frames in frames .

8 Polygons and circular polygons

In this article we have computed conformal maps by solving Laplace problems via least-squares fitting on the boundary by the “lightning” method of lightning , involving poles exponentially clustered outside . Thus we do not make use of explicit analysis of singularities, as is the more familiar approach in numerical conformal mapping. For polygons, the famous method based on analysis of singularities is the Schwarz–Christoffel formula SCbook , which is realized numerically in the widely used and wonderfully robust SC Toolbox by Driscoll toolbox . For circular polygons, there is an analogous procedure based on the Schwarzian differential equation, though it has had less impact numerically bg ; howell .

For both polygons and circular polygons, an in-between option could be employed that would have some advantages. Locally near a corner of a polygon with interior angle for some the conformal map will be an analytic function of (proof by the Schwarz reflection principle). A similar formula applies at a corner of a circular polygon, derivable by a Möbius transformation to a sector bounded by two straight lines. Thus in both of these cases, one can dispense with exponentially clustered poles outside and use a series in terms of explicit singularities instead. At the corner of the polygon, for example, the terms will be By constructing the matrix from columns corresponding to such terms, one can set up a least-squares problem just as before, but making use of a very different set of basis functions, and the convergence will be exponential rather than root-exponential. This would be a numerical method midway between the traditional Schwarz–Christoffel or Schwarzian differential approaches and the fully general one we have focused on in this paper. (A disadvantage is that the inverse map would have to be treated separately, if AAA rational approximation is to be avoided.) We have not pursued this idea, but proofs of concept for L-shaped and circular L-shaped regions, with short Matlab codes, can be found in example .

One might suppose that these observations indicate that there is no need for lightning Laplace solvers in general when the domain is a polygon or circular polygon. However, this is not so, for the Laplace problems that arise in conformal mapping are special ones, with corner singularities determined only by the local geometry lehman . More general Laplace problems often have different boundary data on different sides of a domain, and in such a case, the corner singularities may be complicated even when the geometry is simple lehman2 .

9 Discussion

The algorithms proposed in this paper combine two phases. First, a Laplace problem is solved by least-squares collocation on the boundary. This is essentially the idea of the Method of Fundamental Solutions, which has been pursued by Amano and his coauthors for conformal mapping amano ; amano2 under the name of the charge simulation method. However, our approach differs in the use of poles exponentially clustered near corner singularities to get root-exponential convergence, as proposed in lightning (or in Section 8, the use of explicit singular terms), and in the use of rational functions (dipoles) rather than point charges (monopoles). Second, the resulting map and its inverse are both approximated by rational functions, as proposed in gt .

The idea of solving Laplace problems by fitting Dirichlet data on the boundary is an old one, and early discussions of note are by Walsh walsh , Curtiss curtiss , Collatz (collatz, , sec. V.4.2), and Gaier (gaier, , sec. 4.1), the last of these specifically for conformal mapping. However, these authors formulated systems of equations rather than least-squares problems, making the method fragile and impractical.222Gaier: “Since the configuration of the maximal system of points on the boundary is unknown (except for circles and ellipses), this method is mainly of theoretical interest.” Curtis: “The success of this method is highly sensitive to the correct placement of the points.” Collatz: “The choice of collocation points is a matter of some uncertainty.” As soon as one moves from square matrices and systems of equations to rectangular matrices and least-squares, the distribution of sample points largely ceases to be a problem.

Our methods are based on a minimum of geometrical analysis, and we do not claim they are optimal for any particular problem. For mapping polygons, in particular, the SC Toolbox is tried and true and able to handle quite extreme geometries without difficulty toolbox . Nevertheless, the approach to conformal mapping presented here is very general and usually very fast, enabling typical problems to be solved in on the order of a second of computer time, with subsequent evaluations of both the map and its inverse in microseconds. Experimental codes can be found at lightcode . Extensions to problems such as domains with slits, unbounded domains, or domains bounded by arcs other than straight lines or arcs of circles should be straightforward, though they have not yet been explored.

I have benefited from helpful comments of Toby Driscoll, Abi Gopal, and Yuji Nakatsukasa.


  • (1) B. Adcock and D. Huybrechs, Frames and numerical approximation, SIAM Rev., 61 (2019), pp. 443–473.
  • (2) K. Amano, A charge simulation method for the numerical conformal mapping of interior, exterior and doubly-connected domains, J. Comp. Appl. Math., 53 (1994), pp. 353–370.
  • (3) K. Amano, D. Okano, H. Ogata, and M. Sugihara, Numerical conformal mappings onto the linear slit domain, Japan J. Indust. Appl. Math., 29 (2012), pp. 165–186.
  • (4) S. Axler, Harmonic functions from a complex anlaysis viewpoint, Amer. Math. Monthly, 39 (1986), pp. 246–258.
  • (5) M. Badreddine, T. K. DeLillo, and S. Sahraei, A comparison of some numerical conformal mapping methods for simply and multiply connected domains, Discret. Contin. Dyn. Syst. Ser. B, 24 (2019), pp. 55–82.
  • (6) A. H. Barnett and T. Betcke, Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains, J. Comput. Phys., 227 (2008), pp. 7003–7026.
  • (7) B. Beckermann and A. Townsend, Bounds on the sngular valuces of matrices with displacement structure, SIAM Rev., 61 (2019), pp. 319–344.
  • (8) C. J. Bishop, Conformal welding and Koebe’s theorem, Annals Math., 166 (2007), pp. 613–656.
  • (9) P. Bjørstad and E. Grosse, Conformal mapping of cirulcar arc polygons, SIAM J. Sci. Comp. 8 (1987), pp. 19–32.
  • (10) L. Collatz, The Numerical Treatment of Differential Equations, 3rd ed., Springer, 1960.
  • (11) Computational Methods and Function Theory, special issue on numerical conformal mapping, vol. 11, no. 2 (2012), pp. 375–787.
  • (12) D. Crowdy, The Schwarz–Christoffel mapping to bounded multiply connected polygonal domains, Proc. Roy. Soc. A, 461 (2005), pp. 2653–2678.
  • (13) J. H. Curtiss, Solutions of the Dirichlet problem in the plane by approximation with Faber polynomials, SIAM J. Numer. Anal., 3 (1966), pp. 204–228.
  • (14) H. D. Däppen, Die Schwarz–Christoffel-Abbildung für zweifach zusammenhängende Gebiete mit Anwendungen, diss. ETH Zurich, 1988.
  • (15) T. K. DeLillo, A. R. Elcrat, E. H. Kropf, and J. A. Pfaltzgraff, Efficient calculation of Schwarz–Christoffel transformations for multiply connected domains using Laurent series, Comput. Meth. Funct. Th., 13 (2013), pp. 307–336.
  • (16) T. K. DeLillo, A. R. Elcrat, and J. A. Pfaltzgraff, Schwarz–Christoffel mapping of multiply connected domains, J. d’Anal. Math., 94 (2004), pp. 17–47.
  • (17) T. A. Driscoll, Algorithm 756: A MATLAB toolbox for Schwarz–Christoffel mapping, ACM Trans. Math. Softw., 22 (1996), pp. 168–186; see also
  • (18) T. A. Driscoll, N. Hale, and L. N. Trefethen, eds., Chebfun User’s Guide, Pafnuty Publications, Oxford, 2014; see also
  • (19) T. A. Driscoll and L. N. Trefethen, Schwarz–Christoffel Mapping, Cambridge U. Press, 2002.
  • (20) G. Fairweather and A. Karageorghis, The method of fundamental solutions for elliptic boundary value problems, Adv. Comput. Math., 9 (1998), pp. 69–95.
  • (21) D. Gaier, Konstruktive Methoden der konformen Abbildung, Springer, 1964.
  • (22) D. Gaier, Lectures on Complex Approximation, Birkhäuser, 1987.
  • (23) A. Gopal and L. N. Trefethen, Representation of conformal maps by rational functions, Numer. Math., 142 (2019), pp. 359–382.
  • (24) A. Gopal and L. N. Trefethen, Solving Laplace problems with corner singularities via rational functions, SIAM J. Numer. Anal., 57 (2019), pp. 2074–2094.
  • (25) P. Henrici, Applied and Computational Complex Analysis. III, Wiley, 1974.
  • (26) L. H. Howell, Numerical conformal mapping of circular arc polygons, J. Comp. Appl. Math., 46 (1993), pp. 7–28.
  • (27) C. Hu, Algorithm 785: A software package for computing Schwarz–Christoffel conformal transformation for doubly connected polygonal regions, ACM Trans. Math. Softw., 24 (1998), pp. 317–333.
  • (28) N. Kerzman and E. M. Stein, The Cauchy kernel, the Szegő kernel, and the Riemann mapping function, Math. Annal., 236 (10978), pp. 85–93.
  • (29) N. Kerzman and M. R. Trummer, Numerical conformal mapping via the Szegő kernel, J. Comp. Appl. Math., 14 (1986), pp. 111–123.
  • (30) R. S. Lehman, Development of the mapping function at an analytic corner, Pacific J. Math., 7 (1957), pp. 1437–1449.
  • (31) R. S. Lehman, Developments at an analytic corner of solutions of elliptic partial differential equations, J. Math. Mech., 8 (1959), pp. 727–760.
  • (32) Y. Nakatsukasa, O. Sète, and L. N. Trefethen, The AAA algorithm for rational approximation, SIAM J. Sci. Comp., 40 (2018), pp. A1494–A1522.
  • (33) Y. Nakatsukasa and L. N. Trefethen, An algorithm for real and complex rational minimax approximation, SIAM J. Sci. Comp., submitted.
  • (34) Y. Nakatsukasa and L. N. Trefethen, Vandermonde meets Arnoldi, manuscript in preparation.
  • (35) M. M. S. Nasser, PlgCirMap: A MATLAB toolbox for computing the conformal mapping from polygonal multiply connected domains onto circular domains, arXiv:1911.01787, 2019.
  • (36) D. J. Newman, Rational approximation to , Mich. Math. J., 11 (1964), pp. 11–14.
  • (37) M. Schiffer, Some recent developments in the theory of conformal mapping, appendix to R. Courant, Dirichlet’s Principle, Conformal Mapping, and Minimal Surfaces, Interscience, 1950.
  • (38) H. Stahl, Spurious poles in Padé approximation, J. Comput. Appl. Math., 99 (1998), pp. 511–527.
  • (39) A. Townsend and H. Wilber,

    On the singular values of matrices with large displacement rank,

    Lin. Alg. Applics., 548 (2018), pp. 19–41.
  • (40) L. N. Trefethen, ed., Numerical Conformal Mapping, North-Holland, 1986.
  • (41) L. N. Trefethen, Series solution of Laplace problems, ANZIAM J., 60 (2018), pp. 1–26.
  • (42) L. N. Trefethen, Conformal mapping of L-shaped regions, Chebfun example at, October 2019.
  • (43) L. N. Trefethen, Lightning Laplace and conformal mapping codes laplace.m and conformalR.m,, 2019.
  • (44) J. L. Walsh, Interpolation and Approximation by Rational Functions in the Complex Domain, AMS, 1935.
  • (45) R. Wegmann, Methods for numerical conformal mapping, in Handbook of Complex Analysis: Geometric Function Theory, v. 2, R. Kühnau, ed., Elsevier, 2005, pp. 351–477.