A Comparison of some numerical conformal mapping methods for simply and multiply connected domains

This paper compares some methods for computing conformal maps from simply and multiply connected domains bounded by circles to target domains bounded by smooth curves and curves with corners. We discuss the use of explicit preliminary maps, including the osculation method of Grassmann, to first conformally map the target domain to a more nearly circular domain. The Fourier series method due to Fornberg and its generalization to multiply connected domains are then applied to compute the maps to the nearly circular domains. The final map is represented as a composition of the Fourier/Laurent series with the inverted explicit preliminary maps. A method for systematically removing corners with power maps is also implemented and composed with the Fornberg maps. The use of explict maps has been suggested often in the past, but has rarely been carefully studied, especially for the multiply connected case. Using Fourier series to represent conformal maps from domains bounded by circles to more general domains has certain computational advantages, such as the use of fast methods. However, if the target domain has elongated sections or corners, the mapping problems can suffer from severe ill-conditioning or loss of accuracy. The purpose of this paper is to illustrate some of these practical possibilites and limitations.


1.
Introduction. This paper compares some methods for computing conformal maps from simply and multiply connected domains bounded by circles to target domains bounded by smooth curves and curves with corners. We discuss the use of explicit preliminary maps, including the osculation method of Grassmann [17], to conformally map the target domain to a more nearly circular domain. The Fourier series method due to Fornberg [15] and its generalization to multiply connected domains [3] are then applied to compute the maps to the nearly circular domains. The final map is represented as a composition of the Fourier/Laurent series with the inverted explicit preliminary maps. A novel method for systematically removing corners with power maps is also implemented and composed with the Fornberg maps (which require smooth boundaries) and the level of error that can be expected when using Fourier series to treat domains with corners is illustrated by comparing the results with the Schwarz-Christoffel mapping from SC Toolbox [14] . and the ellipse successively to conformally map the target domain to a more nearly circular domain. The process is a composition of Fornberg maps and inverted Grassmann maps g similar to the simply connected cases in Section 2. The second part covers a method for removing corners in exterior multiply connected domains. The Koebe-like method from Section 2 is used to remove multiple corners from each boundary curve of three rectangles after first inverting in a circle; see Figure 24, below.
In Section 4, we discuss the classical Karman-Trefftz transformation k for removing single trailing-edge corners from an airfoil domain and calculating potential flow. The final map is then k −1 • h, where h is the Fornberg map to the smooth domain. Flow over a single-element and a two-element airfoil are shown, similar to the setup in [21].
Section 5 briefly discusses the bounded, multiply connected case and Section 6 states some conclusions and directions for future research. , where Γ is the inverted ellipse with α = 0.2 (right), produced with composition of 10 Grassmann maps g = g 10 • · · · g 2 • g 1 .
2. Simply connected domains. We recall the Riemann mapping theorem for a simply connected domain.
Theorem 2.1. Let Ω be a simply connected region which is not the whole plane or the Riemann sphere, and let a be a point of Ω. Then there exists in Ω a unique analytic function g satisfying the conditions g(a) = 0 and g (a) > 0, and assuming every value in the unit disk D = {w : |w| < 1} exactly once.
Recall that for conformal maps, g (z) = 0, so that the directions of the curves passing through any point are unchanged. For the proof of the mapping theorem, see, e.g., [1, pp. 325-326]. We will mainly be interested in numerically approximating the inverse of the Riemann map, f = g −1 . (In fact, the composition of the explicit maps can be viewed as an approximation to g.) 2.1. Fornberg's method for the disk. Fornberg [15] proposed a method for computing the conformal map f from the interior of the unit disk to the interior of a smooth, closed curve. We use a slight modification of Fornberg's original formulation which we review briefly here; see [8,13,43,44] for details.
We want to find the conformal map from the interior of unit disk D onto the interior Ω of a Jordan curve Γ. Consider the boundary Γ of Ω to be parametrized by S (e.g., arclength or polar angle), Γ : γ(S), 0 ≤ S ≤ L, γ(0) = γ(L). The normalization imposed on f is f (0) = a ∈ Ω and f (1) = γ(0). Finding f is equivalent of finding the boundary correspondence, S = S(θ) such that f (e iθ ) = γ(S(θ)). If f has the expansion f (e iθ ) = see [24,Sec. 14.3.1]. A Newton-like method is used for numerically solving for the boundary correspondence function S(θ). Let S k (θ) be an approximation to the boundary correspondence function S(θ) at the kth Newton step. We need to find a 2π-periodic correction U k (θ) such that f (e iθ ) = γ(S k (θ) + U k (θ)) are the boundary values at e iθ of the conformal mapping f . Since it is difficult to find the desired correction U k (θ), we compute the correction by linearizing about S k (θ), The condition that the values should have negatively indexed Taylor coefficients leads to a linear problem for the unknown U k . The system is discretized with N -point trigonometric interpolation and results in a symmetric positive definite system, where A is the discretization the identity plus a low rank operator. The system can be solved efficiently with the conjugate gradient method with matrix-vector multiplication using the fft at a cost of O(N log N ). The Newton update is then and near-quadratic convergence is generally observed for a sufficiently close initial guess; see Tables 1 and 5. Usually S 0 (θ) = Lθ/2π will work for domains that aren't too far from circular. Once the boundary correspondence is found, we can easily find the Taylor series, Therefore the Taylor coefficients of f (z) are the Fourier coefficients of the 2πperiodic function, γ(S(θ)), a k := 1 2π 2π 0 γ(S(θ))e −ikθ dθ.
NUMERICAL CONFORMAL MAPPING METHODS 59 2.2. Parametrizing the boundary. Piecewise smooth boundaries may be given by analytic formulas. However, the preliminary maps are constructed by successively conformally transforming a finite set of boundary points. The resulting curve must be parametrized in order to apply Fornberg's method. We fit the points by two periodic cubic splines-one for the x coordinates and one for the y coordinates-parametrized by the chordal arclength between two successive points (x k , y k ), k = 1, . . . , N s along the boundary. Our Matlab code is based on [26]. A smoother fit to the boundary points may be computed with trigonometric interpolation. However, as we show in an example below, this does not to result in much improvement in our computations.

2.3.
Numerical examples using Fornberg's method for analytic boundaries. In this section, we apply the Fornberg map directly and in combination with preliminary Grassmann maps to some explicity known test cases. The calculated errors for the direct Fornberg (Taylor series) map for analytic parametrization of the boundary exhibit spectral accuracy, that is, doubling the number of Fourier points N (or equivalently, the number of Taylor coefficients N/2), squares the error, as expected, since the Taylor coefficients decay geometrically. We consider two examples of test cases for families of curves: the family of inverted ellipses and the family of ellipses; see [7] for other test cases. Each family has a parameter, such as the "thinness" 0 < α ≤ 1 of the domain, that can be changed to make a more poorly conditioned problem requiring larger N to maintain accuracy as α → 0; see [7] for details. Fornberg's method may not converge using the standard initial guess, whereas, Fornberg maps to near circles produced by preliminary Grassmann maps don't require special initial guesses. One could use the preliminary map procedure to generate a good initial guess to Fornberg's direct method for difficult domains, but we have not tried this.
Example 1 (Inverted Ellipse). The boundary curve for the inverted ellipse is 0 < α < 1 is the distance of 0 to the nearest boundary point. The exact map is Note: This map can be derived from the Joukowski map f (z) = z + 1/z which maps exteriors of circles to exteriors of ellipses by inverting, normalizing properly, and rotating. The images of a polar coordinate grid in the interior of the disk under the conformal mapping is displayed in Figure 2. In this case Fornberg's method usually converges with 10 Newton iterations or fewer when the initial guess for the boundary correspondence is given as the standard initial guess, S 0 (θ) := Lθ/2π, where here L = 2π. Quadratic convergence for this analytic boundary parametrization is shown in Table 1. The convergence rate is generally somewhat slower for spline boundaries; see Table 5. The errors in Figure 3 for an analytic parametrization of the boundary exhibit spectral accuracy and are better than the results using preliminary Grassmann maps, Figure 11. However, for α = 0.1 and N = 64, 128, Fornberg's method did not converge using the standard initial guess, whereas, Fornberg maps to near circles produced by preliminary Grassmann maps don't require special initial guesses.
In practice, a boundary curve might be represented by a finite set of points fitted with a periodic cubic spline and spectral accuracy will not be a possibility. Boundary points can be fitted using trigonometric interpolation in an effort to achieve higher accuracy; see, e.g., [37,Fig. 13]. Table 2 gives errors and timings for inverted ellipses, with 10 preliminary Grassmann maps, comparing accuracy and timings for trigonometric interpolation and cubic spline interpolation of the N s nearly circular mapped points. There is only a slight gain in accuracy in some cases using trigonometric interpolation. The points interpolated are not equidistant. Also, the interpolant must be evaluated at the N Fourier points at each Newton iteration of the Fornberg map to the near-circles. For the cubic spline interpolant, this costs O(N ) flops. Since the conjugate gradient iterations for the inner linear systems for Fornberg's method generally converges superlinearly, Fornberg's method costs O(N log N ) using splines (or analytic formulas) for γ(S). (In general, computing the direct map to the domain is slightly faster than using preliminary Grassmann maps with spline interpolation.) However, it costs O(N · N s) flops per Newton step to evaluate the N s-point trigonometric interpolant at the N Fourier points, since the fft cannot be used for the non-equidistant points and it would cost more than O(N ), in any case. The timings in Table 2 for trigonometric interpolation are nearly 10 times greater than for spline interpolation. The small advantage of generating a smoother interpolant with trigonometric interpolation of the very unevenly distributed mapped points thus seems to be lost. It is also necessary to do an additional spline fit to generate a good initial guess of equally distributed points for the Fornberg map to the near circle, since equally distributed points in the trigonometric parameter are not equally distributed on the near circle. If this is not done, the Fornberg iteration will usually not converge. If higher accuracy is needed using the Grassmann maps, it can be achieved by taking large values of N s (and N ) as shown in Figure 12. Table 1. Convergence of successive iteration errors at the Fourier points for the map to an inverted ellipse with α = 0.2 and N = 128. Note that the discrete problem is solved to machine precision but the truncation error is limited by N . Example 2 (Ellipse). See Figure 4. This is an example of an elongated region analyzed in [7]. Maps from the disk to an elongated region have distortions that increase exponentially with the aspect ratio of the region causing very severe ill-condtioning and limiting the applicability of Fourier series methods. It was originally hoped that the use of preliminary maps might circumvent this ill-conditioning, but the large distortions are transferred to the preliminary maps which then amplify any   errors in the Fourier series map to the near circle. Two crescent maps are generally enough to map ellipse domains to very-nearly circular domains. For our calculations, we use the analytic, starlike parametrization of the boundary of the ellipse, γ(S) = ρ(S)e iS , where ρ(S) = α/ 1 − (1 − α 2 ) cos 2 (S) and 0 < α ≤ 1. However, for thin ellipses, using the parametrization γ(S) = cos(S) + iα sin(S) or using a spline fit with, say N s = 2000 knots makes little difference. Note that the conformal map satisfies f (0) = 0, f (±1) = ±1, and f (±i) = ±iα. In [7], we showed that This is an example of the severe ill-conditioning of the conformal mapping problem due to the elongated shape known as the "crowding phenomenon" [35], since images of points at the ends of the ellipse under f −1 are crowded together on the unit circle exponentially in the aspect ratio of the domain as α ↓ 0. For α approaching about 0.2 a large number of Fourier points N are required to achieve any resolution near the ends of the ellipse, if, indeed, the method converges at all; see Figures 15 and 16. For "pinched" domains, such as the inverted ellipses above, the crowding is generally algebraic in 1/α. Note that this ill-conditioning is due to overall shape and not to the presence of corners or high curvature, as discussed in [7], and represents a severe limitation on the use of Fourier series and the unit disk as a computational domain; see, e.g., [11] for an alternative. This ill-conditioning has limited the use of these methods for applications, such as the computation of breaking waves [34] or Rayleigh-Taylor instabilities [35], where highly elongated geometries evolve. In this sense, methods mapping from the physical domain to the disk are much more robust; see [36,38].
Example 3. In practice, there may not be an analytic formula defining the boundary. In Figure 5, the boundary is formed by fitting 10 points in the plane distributed around the origin with a periodic cubic spline using N s = 100 points to parametrize the final curve. The Grassmann maps below can also be applied to this example, but the results look much the same. Note that here the level of error can be no better than the accuracy of the cubic spline used to determine the curve.  The domain Ω to be mapped must be normalized so that it is contained in the unit disk with 0 ∈ Ω and at least two boundary points of modulus 1. The elementary maps gradually expand the domain conformally to fill the disk producing a nearly circular domain. These maps are in the class of osculation maps g i and their composition g k • · · · • g 2 • g 1 converges to the Riemann map at an very slow asymptotic rate of O(1/k). However, they exhibit fast initial convergence as observed by Grassmann and others and analyzed by Henrici [23,24]. In our examples below, k ≤ 10 is generally sufficient to produce a nearly circular domain which can be mapped by Fornberg's method.
We illustrate the behaviour of the Grassmann maps on our two test cases. We will compare the use of inverted Grassmann maps composed with Fornberg's map to the near-circular domain with the computation of the Fornberg map directly to the domain. The first example applies the Grassmann method to map an inverted ellipse to a nearly circular curve after several iterations. The second example is the map to an ellipse. Results for other examples from [7], such as the Cassini oval and the arctanh region yield similar results. Fornberg maps to near circles produced by preliminary Grassmann maps don't require special initial guesses.
We consider two of Grassmann's elementary mappings, the "circle" or "crescent mapping" and the "Koebe mapping". Grassmann used a third elementary map based on the Koebe function, k(z) = z/(1−z) 2 and Porter [40] suggested other ways to accelerate the initial convergence. However, we do not need these alternatives here.
To choose between the two maps at each iteration, we first find the point of modulus A on the boundary of Ω closest to 0. Next, we find the circle K tangent to the boundary of Ω at that point. Ω is rotated so that the point of modulus A goes to −A. If the unit disk containing Ω and the circle K intersect, they form a crescent containing Ω. The crescent can be mapped to unit disk, expanding the point of modulus A to the unit circle.
If K does not intersect the unit disk, then we use the Koebe mapping. Again, Ω is rotated so that the nearest point to the origin maps to −A. Then −A is mapped to the origin by a self-map of the disk, a square root is taken, and a self-map of the disk maps the origin to − √ A. Since 0 < A < 1, we have − √ A < −A, expanding the boundary out toward the unit circle.
2.4.1. Grassmann map calculation. We give a short review of some of the details in [17]. Let z j , j = 1, . . . , N s be points on the boundary Γ of Ω. Let z i0 ∈ Γ be the point on the boundary of Ω nearest to 0. Let where −π < α j < π. The approximate tangent circle K will be the circle through z j , z i0 , z i1 such the Im A j = sin α j is smallest. The calculations of the radius R and the center C of circle K are given by If Im A j < 0, i.e., α j < 0, we delete the interior of K. Otherwise, if Im A j > 0, i.e,. α j > 0, we delete the exterior of K. We multiply the region Ω and circle K by C/ |C| so that the center lies on the x-axis. If K intersects the unit circle in two points we will use the crescent map. Otherwise, we use the second Koebe routine. The intersection points between K and unit circle are given by where Z 1 and Z 1 are the calculated points of intersection.

2.4.2.
Circle or crescent map. Let K be the disk that intersects the boundary of the simply connected curve Ω. The crescent map will be constructed by the composition of the following elementary mappings. Let Z 1 and Z 1 be the intersection points between K and Ω from the previous calculations, and let z j0 be the intersection of circle K and the real axis. The first elementary transformation of the crescent map is the Möbius transformation that will map the crescent to a wedge with angle φ, The second elementary transformation maps to the half plane by raising the first elementary transformation to the power π/φ Finally, the half plane is mapped to the disk by the elementary transformation Therefore, the crescent map will be the composition of the all three elementary transformations, Figure 6 shows the result of a crescent map applied to an ellipse, where the center of the circle tangent to the ellipse at the point nearest the origin lies to the right of the origin. Figure 7 shows the result of a crescent map applied to an inverted ellipse, where the center of the circle tangent to the inverted ellipse at the point nearest the origin lies to the left of the origin.  [1,24]. It moves the boundary to the unit circle more slowly than the crescent map, but it can be applied to any boundary. Here the tangent circle K is entirely in the unit disk where z i0 , the point closest to 0, is rotated to −A = −|z i0 |. R is the radius and C is the center of the circle tangent to the boundary at −A. This means C < 0 and R + |C| < 1. We map the unit disk slit along negative x-axis to −P = −(C + γR) where 0 < γ < 1 and −1 < −P < 0 since 0 ∈ Ω. Here we use γ = 0.9, so that square roots aren't taken on the boundary. The first elementary mapping is to map the point −P to the center 0 by a self-map of the unit disk, z − − → w = z + P 1 + P z . The next elementary transformation is the square root that maps the slit to the half disk, w − − → √ w, and the last transformation to the desired shape is The Koebe map will be the composition of the all three elementary transformations, Note that −P → 0 → − √ P < −P and so that the boundary near −A ≈ −P will be moved toward the unit circle. Figure 8 illustrates the steps in the construction for an inverted ellipse where the tangent circle does not intersect the unit circle and the crescent map cannot be used.

2.4.4.
Inverse crescent map. The inverse of the crescent mapping uses the following conformal mapping transformations.
Inverse Koebe map. The inverse of the Koebe mapping uses the following conformal mapping transformations.
Example 4 (Grassmann iterations). Some iterations of the Grassmann maps applied to an inverted ellipse, Example 1 with "thinness" α = 0.2, are shown in Figures 9 and the rate of convergence of the boundary to the unit circle for various α is illustrated in Table 3. Note that the initial convergence rate is fast, as discussed in [23,24,40]. The resulting distribution of N s = 400 points on the near-circle after igm := the number of Grassmann maps = 10 is shown in Figure 10. The boundary of the nearly circular domain is parametrized by fitting a periodic cubic spline parametrized by chordal arc length, as descibed above, through these points. To increase the overall accuracy more spline points N s can be mapped at a relatively  Figure 8. Steps for one Koebe map applied to an inverted ellipse with α = 0.48. Note that the circle tangent to the boundary at the point nearest the origin does not intersect the unit circle, so that the crescent map cannot be applied. The square root is taken at the origin slightly off of the boundary curve, so that the curve remains smooth. cheap cost of order N s× the number of Grassmann iterations. Examples of other curve families, such as those discussed in [7], have been computed in [2,4] with similar results.  Figure 9. Iterations 1, 3, and 6 of Grassmann maps applied to an inverted ellipse with α = 0.2; see Table 3.   2.5. Removing corners. In this section, we introduce a method for systematically removing corners. The method is based on the Koebe map above for general curves. An example of a square is used to calculate the error in removing its corners with our  Koebe-like method compared to the Schwarz-Christoffel map using SC Toolbox [14], which is highly accurate. See [9,10,18,29,30,31,44] for related papers. The domain Ω is normalized inside the unit disk such that 0 ∈ Ω and such that none of the corners lie on the unit circle. We apply the following procedure, similar to the  Koebe map above, successively to each corner. First, we rotate the corner z k with |z k | = A < 1 to −A. Next we apply the following sequence of maps, where β k π is the corner angle, One step of this procedure is shown in Figure 17.  Figure 19. The error is computed in Table 4 by comparison with the map produced by SC Toolbox [14]. Even with N = 1024 Fourier points, we only get about 3 digit accuracy. This is consistent with our other experience computing maps to domains with corners with methods based on Fourier series in [9,10]. Typically one may see small oscillations of the boundary near the corners. (A careful study of the behaviour of Theodorsen's method, a Fourier series method, combined with techniques for smoothing the series, such as Lanczos smoothing, for domains with corners is given in [18].) A map to a polygon with 5 corners, including a re-entrant corner, is shown in Figure 20. Also, a map to a crescent is shown in Figure 21. In this case the exact map is known and the accuracy is only about 2 digits. This is a simple example of a curvilinear domain with corners. SC Toolbox is restricted to polygonal domains. However, a "generalization" of the Schwarz-Christoffel transformation to exterior curvilinear polygons was considered in [10], but did not give high accuracy.
To treat curvilinear domains with high accuracy, methods which map to the disk and can refine the mesh in the target domain near the corners should be used; see, e.g., [36] and also [39,Remark 1.5.3] and [19,20] for alternative methods which can treat corners and even cusps using singular functions or adaptive methods.  Figure 17. The steps for smoothing a corner: The domain is placed slightly inside of the unit disk and each corner is successively rotated to the negative real axis and mapped to the origin with a Möbius map of the disk taking the corner to the origin. The corner is smoothed with a power map, and the origin is mapped back to the negative real axis.       3.1. Extension of Fornberg's method to exterior multiply connected domains. The details of this method are presented in [3]. We review them briefly here.
Boundaries of disks D k are the circles C k : c k (θ) := z k +ρ k e iθ and C = C 1 +· · ·+C m . The target boundary of Ω is Γ = Γ 1 + · · · + Γ m , where the Γ k are (smooth) curves parametrized by S, e.g., arclength, Γ k : γ k (S). f extends smoothly to the boundary f (C k ) = Γ k . To compute the map, we must find the boundary correspondences S = S k (θ) and conformal moduli z k , ρ k such that where f is analytic and f (z) = z + O(1/z), z ≈ ∞. Therefore, As in the simply connected case, we will linearize the problem about the current guesses for the boundary correspondences and centers and radii and use a Newtonlike iteration. For an initial guess S k (θ) and 2π periodic correction U k (θ), we use the following linearization for U k (θ) For an initial guess of z k and ρ k with corrections δz k and δρ k , The above approximations gives Requiring the functions to be boundary values of a function analytic in D gives a linear system for the unknown Newton updates U k , δz k , δρ k denoted by U , of the form, AU = b, where A is a symmetric positive definite matrix of the form identity plus a low rank operator. The system can be solved efficiently by the conjugate gradient method, as in the simply connected case, except that the matrix-vector multiplication now costs O((mN ) 2 ) operations instead of O (N log N ).
The i + 1st Newton updates are then, Example 6 (Preliminary Grassmann maps for multiply connected domains). In this section, we apply Grassmann maps to multiply connected exterior domains and compose them with the Fornberg maps. Figure 22 shows a domain of connectivity m = 4 mapped with Grassmann maps to produce the near circles. The domains are successively inverted and the Grassmann maps for the interior are used to produce the near circles. The extension of Fornberg's method to the exterior multiply connected case [3] is used to compute the map from the exterior of m = 4 circles to the exterior of the near circles and then the inverse Grassmann maps map to the original domain. The differences in the successive iterates of the Fornberg method are shown in Table 5 and indicate that the discrete problems is solved accurately. Fornberg's method to map directly to this domain did not converge using the standard initial guess. For the domain in Figure 23, both the direct map to the domain and the map computed with Grassmann maps converged. Table 5. Convergence of successive Newton-iteration errors S i+1 − S i ∞ at the Fourier points for the Fourier series map to the the near-circular region of connectivity m = 4 in Figure 22 for N = 128, 256. Note that the convergence rate is nearly independent of N . ( 3.1e-14 5.6e-14 10 1.8e-14 2.7e-14 Table 6. Some typical sample timings for map to regions in Fig  In Figure 24 the Fornberg map is composed with successive applications of our corner smoothing method to produce the map to the exterior of m = 3 rectangles. Table 7 shows the convergence of the successive iterates for near circular domains resulting from domains with corners.
4. Karman-Trefftz transformations and potential flow. A popular method for removing a single corner on the boundary curve of a domain exterior to the curve and containing ∞ is the Karman-Trefftz transformation. This has been used often for computing potential flow over the exterior of a multi-element airfoil [21]. If the exterior angle is then βπ, the Karman-Trefftz transformation is given by where z 1 , the corner at the trailing edge, maps to ζ 1 and z 2 , a point chosen interior to the curve near the leading edge, maps to ζ 2 . Then and its inverse is The application of k(z) to the airfoil usually results in a nearly circular set of point which we fit with our periodic cubic spline routine parametrized by chordal  arclength, as in the application of the osculation maps. For the multiply connected case, the Karman-Trefftz transformation can be applied successively to the images of the m airfoils to produce a map k = k m • · · · • k 2 • k 1 from the domain Ω bounded by airfoils to the domain k(Ω) bounded by nearly circular, smooth curves. Note that each of these maps k i must be applied to all of the curves at each step. Fornberg's method for exterior multiply connected domains [3] can then be used to compute a Laurent series map h from a conformally equivalent domain D exterior to m disks to k(Ω). The Karman-Trefftz maps k i can be explicitly inverted. The final conformal map f from the circle domain D to the domain Ω exterior to the m airfoils can be represented as a composition,  Figure 26 for a single Joukowski airfoil and Figure 28 for m = 2 cosine curves. An X marks the z 1 , z 2 , ζ 1 , ζ 2 for each (transformed) domain. We compute the streamlines about the airfoils by adding circulation to streaming flow in order to satisfy the Kutta-Joukowski condition at the trailing edge, as in [21,45]. A more complete discussion of this procedure will be given in future work. Example 7. The first example is the map to the exterior of a Joukowski airfoil Figure 26. This is a well-known case where the map from the exterior of a circle to the exterior of an airfoil with given trailing edge angle is known explicitly; see [33]. Specifically, the map used here to an airfoil with interior trailing edge angle 0.15π radians is given by the sequence of maps, Figure 25 illustrates the case where N s = 101 points, equally spaced in θ, are distributed about the airfoil and mapped with k(z) to a smooth domain. We deliberately do not just invert the Joukowski map to the circle in order to have an exact test case. The errors between the exact Joukowski map and k −1 •h for increasing values of N and N s are given in Figure 27. Note that the error can be made as small as we please, since the preimage of the one vertex at the trailing edge can always be made to be a Fourier point. This is not possible, in general, for connectivity greater than one. Example 8. The second example is the domain exterior to two cosine airfoils. The boundaries are given by with a trailing edge interior angle π K and K = 4. We apply this to a distribution of 400 points along the cosine curve with angle βπ = 2π − π/K, and 1/β = K 2K−1 . Table 7 lists the maximum error between the successive iterations of the boundary correspondence functions S i 1 and S i 2 , illustrating the near quadratic convergence of Fornberg's method. (Recall that this does not give the discretization error which depends on N , but only the accuracy of the solution to the discrete problem. This rate is independent of N .) In, e.g., [21] and earlier work, linearly convergent Fourier series methods, such as those of Theordorsen or Timman/James [27] for simply connected domains applied successively to each boundary, were used for the maps to the smoothed domains; see [24] for a presentation of these methods and [9] for comparisons with the methods of Fornberg and Wegmann for simply connected regions. In [9], a discretization error of 10 −4 was achieved with this approach for the simply connected case with N = 256. 5. The bounded, multiply connected case. A version of our Fourier series method for mapping to bounded multiply connected regions with smooth boundaries has been developed in [12]. We have not done much testing of this method with explicit preliminary maps, but we expect the results to be much the same as the unbounded case, above. Figure 29 shows the map to the interior of a diamond with elliptical holes, where the corners have been smoothed using our Koebe-like method. The straight sides are well represented, but it is again difficult to resolve the corners very sharply.  6.4e-02 5.4e-10 1.1e-10 7 6.9e-03 5.4e-10 3.3e-13 8 1.2e-04 5.4e-10 1.7e-15 9 1.9e-06 5.4e-10 4.1e-15 10 4.0e-08 5.4e-10 9.7e-15 11 6.2e-10 12 1.3e-11 13 1.8e-13 14 1.6e-14 6. Concluding remarks. The results of our computations show that preliminary explicit maps are marginally useful and lead to somewhat more robust methods. An accurate initial guess is easily given for the near-circular regions in both the simply and the multiply connected cases. However, the coding becomes somewhat complicated, especially in the multiply connected case, and accuracy can be lost due to the necessity of parametrizing the nearly circular or smoothed boundaries with curves of lower regularity, such as cubic splines. If the original curves are given by spline fits to a finite number of points defining the boundary, the accuracy will be limited from the start. For domains with corners only 3 or 4 digit accuracy can be expected when composing corner smoothing with series methods, except in special cases such as a single airfoil. In future work we plan to compare our results with Figure 29. The map from the interior of the unit disk (left) with two circular holes to the interior of a diamond (right) with two elliptical holes by smoothing corners on the diamond and using N = 256 Fourier points for the map to the smooth region (middle).
linearly convergent projection methods described in [41,42,44] which can be applied directly to regions with corners. Also, more comparison with Wegmann's Newtonlike methods [44] based on solving Riemann-Hilbert problems would be useful and some comparsions were given in [3]. However, we do not expect any significant advantages overall (except that the Newton-like methods are faster), since these methods are all based on Fourier series. Also, we plan more complete experiments on computations of potential over airfoils.