Schwarz–Christoffel Mapping in the Computer Era

Lloyd N. Trefethen and Tobin A. Driscoll

Abstract. Thanks to powerful algorithms and computers, Schwarz– Christoffel mapping is a practical reality. With the ability to compute have come new mathematical ideas. The state of the art is surveyed. 1991 Mathematics Subject Classification: 30C30, 31A05 Keywords and Phrases: conformal mapping, Schwarz–Christoffel formula

1. Introduction. In the past twenty years, because of new algorithms and new computers, Schwarz–Christoffel conformal mapping of polygons has matured to a technology that can be used at the touch of a button. Many authors have contributed to this progress, including Däppen, Davis, Dias, Elcrat, Floryan, Henrici, Hoekstra, Howell, Hu, Reppe, Zemach, and ourselves. The principal SC software tools are the Fortran package SCPACK [15] and its more capable Matlab successor, the Schwarz–Christoffel Toolbox [3]. It is now a routine matter to compute an SC map involving a dozen vertices to ten digits of accuracy in a few seconds on a workstation.

With the power to compute has come the ability to explore. The obvious kind of problem—“here is a polygon; map it onto a disk or a half-plane”—is rarely what one encounters in practice. Instead it is variations on the idea of SC mapping that arise. In this article, after briefly mentioning the algorithmic developments that have made SC mapping possible, we describe four of these variations: oblique derivative problems on polygons; ideal free-streamline flows; the CRDT (cross-ratio Delaunay triangulation) algorithm; and Green’s functions for symmetric multiply connected domains.

2. Numerical Algorithms. Independently around 1869, Schwarz and Christoffel derived their famous formula,

f(z) = A + B ∫_{0}^{z} ∏_{k=1}^{n} (1 − ζ/z_{k})^{−β_{k}} dζ.  (1)

They proved that any conformal map f(z) of the unit disk or the upper-half plane onto a polygon P with n vertices can be written in the form (1) for some constants A, B, {z_{k}}, and {β_{k}} [9,10]. The exponents {β_{k}} are determined by the angles at the vertices of P (exterior turning angles divided by π), but the other parameters are unknown. Translation, scaling, and rotation are accomplished by A and B, and the prevertices {z_{k}}, unknown a priori, are the preimages on the unit circle or the real axis of the vertices of P.

Three computational hurdles arise in implementing (1): finding the unknown parameters, evaluating the integrals, and computing the inverse map. Analytically, one can do very little. The hurdles must be crossed numerically, and the history here is long, in part because this is a topic that every engineer has heard of and may be tempted to solve from scratch. The list of names above includes only some of the more prominent contributors. Two of the earliest contributions in the list were those of Davis [1] and Reppe [13].

We will not give details, but mention a few of the algorithmic ideas that are the basis of SCPACK and the SC Toolbox. Schwarz–Christoffel integrals can be evaluated by an automatic process of compound Gauss–Jacobi quadrature: Gauss– Jacobi formulas handle the singularities at endpoints, and adaptive subdivision of intervals combats the phenomenon of exponential distortions that is universal in conformal mapping (see §5). The unknown parameters can be found by solving a system of nonlinear equations that assert that the side lengths of P are correct. By a change of variables, the ordering conditions among the prevertices can be eliminated to make this system of equations unconstrained, whereupon it can be treated by quasi-Newton iteration or by more specialized techniques. Finally, f^{−1} can be evaluated by a Newton iteration (the derivative f′ is just the integrand of (1)), made robust as necessary by the generation of initial guesses via numerical solution of an ordinary differential equation.

We urge readers to download a copy of the SC Toolbox [3] and give all this a try. Begin by typing scgui to explore the graphical user interface, but remember that everything can also be done by inline Matlab commands. For example,

      plot(diskmap(polygon([2 2+i 1+i 1+2i 2i 0])))

generates a plot of a conformal map of the unit disk onto an L-shaped polygon in five seconds on the first author’s workstation. Changing diskmap to extermap maps the exterior of the same polygon. The command

      plot(rectmap(polygon([2 2+i 1+i 1+2i 2i 0]),[1 4 5 6]))

maps the interior onto a rectangle, whose length-to-width ratio is necessarily the conformal modulus of the L-shaped region, \sqrt{3} . Similar SC Toolbox commands construct maps from a half-plane or an infinite strip and onto generalized polygons with slits, vertices at infinity, or overlapping regions.

3. Oblique Derivative Problems on Polygons. The following problem, considered by Trefethen and Williams [18], arises in queuing theory and in the study of the classical Hall effect in electronics. On each side Γ_{j} of a polygon P, an angle θ_{j} is specified. We seek a non-constant harmonic function u in P, i.e., a solution to Laplace’s equation Δu = 0, that satisfies the condition du/ds = 0 along the direction at angle θ_{j} to the boundary.

Such problems can be solved by conformal mapping as follows. Suppose a function f is found that maps P conformally onto a second polygon Q, whose side lengths are unspecified, with the property that f(Γ_{j}) is oriented at the angle θ_{j} from the vertical. Then all the boundary directions in Q line up vertically, and therefore u(z) = Ref(z) is the function required. Theorems 1–3 of [18] establish that this procedure generates all solutions to the oblique derivative problem.

Figure 1 shows an example of an oblique derivative problem solved in this manner. Though only pictures are presented, it is an easy matter to extract numbers from such a computation accurate to ten or more digits.

P  f  Q

Figure 1. Solution of an oblique derivative problem by conformal mapping. On the L-shaped problem domain P, a bounded harmonic function u is sought with du/ds = 0 at angle π/4 on the left edge, π/3 on the right edge, and π/2 (the usual Neumann condition) elsewhere. The required function is u(z) = Ref(z), where f(z) is a conformal map of P to a trapezoid Q with sides oriented at the prescribed angles from the vertical, and all solutions are of this form, differing only in shift and scale. The preimages of vertical lines in the trapezoid are curves u(z) = const. The two dots show the conformal images of the two vertices in the L that map into degenerate vertices of the trapezoid.

4. Ideal Free-Streamline Flows. A longstanding topic of fluid mechanics is the study of jets, wakes, and cavities, all of which may involve a surface, or in two dimensions a line, across which the flow properties change. In 1868, Helmholtz and Kirchhoff introduced the theory of free streamlines for such problems and proposed the use of complex analysis to find solutions in 2D. Other early contributors were Planck, Joukowski, Réthy, Levi-Civita, Greenhill, and von Mises, and later generations saw major extensions and survey publications by Birkhoff and Zarantonello, Gilbarg, Gurevich, Monakhov, and Wu, among others [12,20].

The classical theory of 2D free-streamline flows has two limitations. The first is that in many cases it omits important aspects of the physics. The second is that the flows in question can be computed analytically for only the simplest geometries. Here, however, it turns out that by a modification of the SC idea, effective algorithms can be devised that are exactly analogous to those used for SC mapping. If one is careful about the physics, as in various papers over the years by J. Keller and by Vanden-Broeck, among others, the results can reveal a great deal about certain flows.

One version of a classical 2D free-streamline flow problem goes like this. A semi-infinite inviscid, incompressible fluid flows in the absence of gravity above a polygonal wall that ends at a point, after which the fluid continues on into free space. Beyond the separation point, the flow is bounded by a curved free streamline on which the boundary condition (because of Bernoulli’s equation) is that the magnitude of the velocity must be constant—say, 1.

The solution can be found as follows. Let z be the spatial variable; the boundary of the problem domain in the z-plane is partly unknown, because of the free streamline. Let w be the velocity potential, inhabiting the upper halfplane, with the origin corresponding to the point of separation. Let ζ = dw/dz be the hodograph or conjugate-velocity variable. The problem is to find an analytic function ζ(w) such that argζ(w) takes prescribed piecewise constant values along (−∞, 0] and has constant modulus |ζ(w)| = 1 along [0, ∞).

−β_{k}π −β_{k}π ζ_{k} f_{k} z_{k} 0 w_{k}  0  1 0 (a) (b)

Figure 2. Comparison of the ideas underlying Schwarz–Christoffel and freestreamline mapping. (a) For an SC map f(z), the derivative f′(z) has piecewise constant argument for z ∈ (−∞, ∞), so it can be written as a product of elementary maps of the upper half-plane onto infinite wedges. (b) For a free-streamline map, ζ′(w) has piecewise constant argument for w ∈ (−∞, 0] and constant modulus for w ∈ [0, ∞), so it can be written as a product of maps of the upper half-plane onto a circular-arc wedges.

The hodograph (ζ) domain is bounded by straight segments and a circular arc, and the classical approach would be to take the logarithm to reduce it to a polygon, which could then be mapped by the SC formula. For all but very simple geometries, however, this method is unworkable because of unknown ordering of prevertices, for the topology of the polygon is not known in advance. Instead, general polygonal boundaries can be handled by the method suggested in Fig. 2, developed by Monakhov [12] and Elcrat and Trefethen [6]. The key idea is to employ a modification of the SC integral (1) in which each term in the product in the integrand is an elementary map onto a bounded circular-arc wedge rather than an unbounded wedge. A numerical example is presented in Figure 3; for more, see [2] and [6].

The method just described for free-streamline flows is a general one, permitting the ready computation, with slight variations, of flows in a wide variety of geometries. Once such a general tool is in hand it is an easy matter to reproduce computations from the past, including those of Kirchhoff (flat plate, 1869),

Figure 3. Example of a free-streamline flow computed by the method of Fig. 2. The free streamline is the curve that separates from the tip of the solid boundary. The equal spacing of the curves on the right reveals that the constantspeed condition has been satisfied.

Rayleigh (inclined plate, 1876), Bobyleff (symmetric wedge, 1881), Michell (slot, 1890), von Mises (funnel, 1917), Chaplygin and Lavrentiev (plate with separation, 1933), Keller (“teapot effect”, 1957), Lin (asymmetrical wedge, 1961), Wu and Wang (symmetric 4-piece wedge, 1964), and Elcrat (plate with spoiler, 1982).

5. The CRDT (Cross-Ratio Delaunay Triangulation) Algorithm. The methods described above are based on what might be called standard SC mapping technology, in which the realization of (1) is achieved by standard “best practice” methods of numerical analysis. It has been recognized since around 1980, however, that such methods fail for highly elongated regions. Suppose, for illustration, that the unit disk is mapped onto a rectangle of aspect ratio L, with 0 mapping to the center of the rectangle. Then the prevertices along the unit circle lie in two pairs separated by intervals that shrink in proportion to exp(−πL/2). For L = 30, for example, adjacent prevertices are separated by only about 10^{−21}. Thus conformal maps are subject to exponential distortions, a phenomenon known as crowding, and in floating point arithmetic, the result is that highly elongated regions cannot be treated by standard methods.

One solution, due to Howell and Trefethen [11], is to dispense with the disk or half-plane and map directly onto a highly elongated domain such as an infinite strip or a long rectangle. Both options are included in the SC Toolbox. For many problems arising in practice, this solves the crowding problem by weakening the effect from exponential to algebraic.

Domains that are elongated in multiple directions, however, require more radically new approaches, and one such, also included in the SC Toolobx, has recently been developed by Driscoll and Vavasis [4]. The idea behind the CRDT algorithm is that no matter how distorted a conformal map may be globally, any portion of it can be made locally well-behaved by some Möbius transformation (az + b)/(cz + d). By composing Möbius transformations, it ought to be possible to represent arbitrarily great distortions by compositions of well-behaved maps.

The CRDT algorithm begins by constructing a Delaunay triangulation of the target polygon P, generally after introducing extra degenerate vertices so that the triangles will be not too slender. The task then is to find the correspondence function between prevertices on the unit circle and vertices on the boundary of P. The Möbius idea is used by formulating the correspondence condition not globally but four vertices at a time, corresponding to two adjacent triangles. A convenient, Möbius-invariant description of the unknown prevertices is furnished by the crossratios of these 4-tuples, defined by

ρ(z_{1}, z_{2}, z_{3}, z_{4}) = \frac{(z_{4} − z_{1})(z_{2} − z_{3})}{(z_{3} − z_{4})(z_{1} − z_{2})},  (2)

which is negative and real when z_{1}, …, z_{4} lie on counterclockwise order on the unit circle. The CRDT algorithm formulates a system of n − 3 nonlinear equations in which the independent variables are the negatives of the logarithms of the crossratios of 4-tuples of prevertices on the unit circle, and the dependent variables are the deviations from their correct values of the absolute values of the logarithms of the cross-ratios of 4-tuples of vertices on P. This system of equations is observed to be very well behaved and easily solved by iteration.

The algorithm sounds complicated, and its elements of computational geometry certainly give it a flavor different from other algorithms in this field, but it has proved remarkably effective. It makes possible the mapping of regions that would have been regarded as impossible a few years ago, except in multiple precision arithmetic. Figure 4 gives an example.

Figure 4. Conformal map of “Emma’s maze” onto a rectangle of aspect ratio 18.2, computed by the CRDT algorithm. The solid curves are the conformal images of straight lines in the rectangle; the dotted curves are the same, but correspond to lines in the rectangle exponentially close to the sides (10^{−2}, 10^{−4}, 10^{−6}, …, 10^{−16}). Because of exponential distortions, the numerical computation of maps like this is far beyond the capabilities of ordinary algorithms.

6. Green’s Functions for Multiply Connected Domains. Our final SC variation represents recent work by Embree and Trefethen [6] based on ideas going back at least to Widom [19]. It is a longstanding dream to generalize SC maps to multiply connected domains, but except in the doubly-connected case treated in successive works by Henrici, Däppen, and Hu, no general results of much practical value are available in this line. Consider, however, the special case of a region P in the complex plane consisting of polygons P_{1}, …, P_{K} that are symmetrically located along the real axis, illustrated in Fig. 5 for an example with K = 3. (In an important special case, the polygons degenerate to intervals along the real axis.) To be specific, suppose we seek the Green’s function for P, the function u(z) harmonic outside P with boundary values u(z) = 0 on the boundary of P and u(z) ∼ log |z| as z → ∞.

Such a Green’s function can be computed by SC mapping. The crucial observation is that the upper half of the problem domain, with segments of the real axis inserted as necessary to provide a complete boundary, is simply-connected. Let g(z) be a conformal map of this region onto a semi-infinite strip with vertices πi, 0, and ∞ and K − 1 slits of indeterminate height and length along the complex interval [0, πi] (Fig. 5(b)). The semi-infinite segments of the real axis map to the sides [0, ∞) and [πi, ∞) of the strip, and the bounded segments between the polygons P_{j} map to the horizontal slits; the boundaries of the polygons themselves map into segments of [0, πi]. The Green’s function required is now given by

u(z) = Reg(z).  (3)

A second conformal transformation may provide further insight. If f(z) = exp(g(z)), then f maps the upper half of the problem domain onto the exterior of the unit disk with protruding spikes (Fig. 5(c)). In terms of this new map, the Green’s function is given by

u(z) = log |f(z)|.  (4)

This Green’s function algorithm has several special features from the SC mapping point of view, both of which arise in certain other applications too. One is that the positions of the slits (spikes) in the semi-infinite strip (exterior of the disk) are unknown a priori, and must be determined as part of the mapping process. Doing so is easy, since this part of the parameter problem enters linearly, as is often the case with slits in SC mapping; but we can view this aspect of the calculation as prototypical of more complicated generalized parameter problems that arise in various applications such as inverse problems [16].

The second special feature of this SC computation is that although the Green’s function u(z) is single-valued, the conformal maps involved in getting it, if viewed in the large, are multiple-valued. Specifically, consider the function f(z). A priori, it maps the upper half of Fig. 5(a) to the upper half of Fig. 5(c), and reflection in the semi-infinite segments of the real axis completes this to a single-valued map of all of (a) to all of (c). Reflection in the bounded segments of the real axis between the polygons that correspond to the spikes in Fig. 5(c), however, is equally justified mathematically. After one such reflection, further reflections become possible, and (a)

(b)  0 πi  g(z)

(c)  exp

Figure 5. Computation of the Green’s function u(z) for a multiply connected domain. (a) Final result, showing four level curves u(z) = const.; the innermost is the critical one at which two connected components first touch. (b) To compute u, the portion of the problem domain in the upper half-plane is transplanted by an SC map g(z) to a semi-infinite slit strip; the upper halves of the level curves of (a) are the preimages of vertical lines in the strip (dashed). (c) The exponential function transplants the slit strip to the upper half of the exterior of a disk with spikes, and the level curves of (a) become preimages of concentric circles. so on, and the ultimate result is that f(z) maps the multiply connected region onto a Riemann surface with (in general) infinitely many sheets [19]. Fortunately, these complications do not matter for the application of computing the Green’s function.

Green’s functions for multiply connected regions have applications to problems of polynomial approximation. When the polygons reduce to intervals we have a problem in digital filtering [14], and the general case relates to beautiful theorems of Fuchs in the 1970s [8]. Suppose in Fig. 5, for example, that three distinct entire functions f_{1}, f_{2}, f_{3} are defined and polynomials p_{n}(z) of increasing degrees n are sought that approximate f_{j} on P_{j} in the uniform norm. Then apart from a small algebraic factor, the approximation errors decrease at the rate exp(−βn), where β is the length of the shortest slit in Fig. 5(b). The heights of the slits give asymptotic information about the proportions of interpolation points on each set P_{j}, with the images of roots of unity in Fig. 5(c) on the boundaries of P_{j} providing near-optimal interpolation points, and the optimal approximants have the property that they converge precisely inside the critical level curve of Fig. 5(a).

7. Further Variants and Applications. We have touched upon only a few developments in Schwarz–Christoffel mapping in the computer era. Among the variants not mentioned are periodic domains, fractals, circular polygons, curved boundaries, gearlike domains, and polygonal Riemann surfaces. Among the applications not mentioned are Faber polynomials, matrix iterations, the KdV equation, electrical resistances and capacitances, magnetostatics, and vortex methods in fluid mechanics. References for some of these problems can be found in [17].

We would like to put in writing our view of the best applications of Schwarz– Christoffel transformations. Many people have the idea that SC methods may be useful for general geometric purposes such as grid generation for finite differencing, perhaps even for domains with curved boundaries approximated by polygons. Our opinion is that whereas such applications are of course sometimes effective, the real excitement of SC mapping lies elsewhere. What is special about the SC formula is that it solves a certain precisely defined problem exactly, delivering a semi-analytic solution dependent only on a small number of parameters. We favor applications where this semi-analytic solution solves the problem of ultimate interest exactly, for in such circumstances, SC methods far outperform general-purpose tools such as adaptive finite elements.

References

[1] R. T. Davis, Numerical methods for coordinate generation based on Schwarz– Christoffel transformations, Proc. 4th AIAA Comp. Fluid Dyn. Conf., 1979.

[2] F. Dias, A. R. Elcrat and L. N. Trefethen, Ideal jet flow in two dimensions, J. Fluid Mech. 185 (1987), 275–288.

[3] T. A. Driscoll, Algorithm 765: A MATLAB toolbox for Schwarz–Christoffel mapping, ACM Trans. Math. Softw. 22 (1996), 168–186. Software available at

http://amath.colorado.edu/appm/faculty/tad/research/sc.html.

[4] T. A. Driscoll and S. A. Vavasis, Numerical conformal mapping using crossratios and Delaunay triangulation, SIAM J. Sci. Comput., to appear.

[5] A. R. Elcrat and L. N. Trefethen, Classical free-streamline flow over a polygonal obstacle, J. Comput. Appl. Math. 14 (1986), 256–265.

[6] M. Embree and L. N. Trefethen, Green’s functions for multiply connected domains via Schwarz–Christoffel mapping, manuscript in preparation.

[7] J. M. Floryan and C. Zemach, Schwarz–Christoffel mappings: A general approach, J. Comput. Phys. 72 (1987), 347–371.

[8] W. H. J. Fuchs, On the degree of Chebyshev approximation on sets with several components, Izv. Akad. Nauk. Armyan. SSR 13 (1978), 396–404.

[9] D. Gaier, Konstruktive Methoden der Konformen Abbildung, Springer, Berlin, 1964.

[10] P. Henrici, Applied and Computational Complex Analysis, v. 3, Wiley, New York, 1986.

[11] L. H. Howell and L. N. Trefethen, A modified Schwarz–Christoffel transformation for elongated regions, SIAM J. Sci. Stat. Comput. 11 (1990), 928–949.

[12] V. N. Monakhov, Boundary-Value Problems with Free Boundaries for Elliptic Systems of Equations, Trans. of Math. Monographs, vol. 57, Amer. Math. Soc., 1983.

[13] K. Reppe, Berechnung von Magnetfeldern mit Hilfe der konformen Abbildung durch numerische Integration der Abbildungsfunktion von Schwarz–Christoffel, Siemen Forsch. u. Entwick. Ber. 8 (1979), 190–195.

[14] J. Shen, G. Strang, and A. Wathen, The potential theory of several intervals and its applications, manuscript in preparation.

[15] L. N. Trefethen, Numerical computation of the Schwarz–Christoffel transformation, SIAM J. Sci. Stat. Comput. 1 (1980), 82–102.

[16] L. N. Trefethen, Analysis and design of polygonal resistors by conformal mapping, J. Appl. Math. Phys. 335 (1984), 692–704.

[17] L. N. Trefethen, Schwarz–Christoffel mapping in the 1980s, TR 93-1381, Dept. of Comp. Sci., Cornell U., 1993.

[18] L. N. Trefethen and R. J. Williams, Conformal mapping solution of Laplace’s equation on a polygon with oblique derivative boundary conditions, J. Comp. Appl. Math. 14 (1986), 227–249.

[19] H. Widom, Extremal polynomials associated with a system of curves in the complex plane, Adv. Math. 3 (1969), 127–232.

[20] T. Y. Wu, Cavity and wake flows, Annual Reviews in Fluid Mech. (1972), 243–84.

Lloyd N. Trefethen

Oxford U. Computing Laboratory Wolfson Building, Parks Road Oxford OX1 3QD, UK

     LNT@comlab.ox.ac.uk

Tobin A. Driscoll

Dept. of Applied Mathematics University of Colorado Boulder, CO 80309-0526, USA

                                      driscoll@na-net.ornl.gov