Michael D. O'Neill

In this article we  reproduce the 1979 experiments of E. Grassmann on approximation of the Riemann mapping by compositions of elementary  conformal  mappings. The experiments were somewhat surprising in that they showed a fast initial speed of the algorithm. Ostrowski had shown  that the rate of   convergence for the method would, in general, be asymptotically slow. The function given here for implementing Grassmann's version of the algorithm provides a nice classroom demonstration of a constructive proof of the Riemann mapping theorem. Such a proof requires the Schwarz lemma, Harnack's inequality and the argument principle but not the theory of normal families.

1. Introduction

The Riemann Mapping Theorem

The Riemann mapping theorem states that any simply connected  region ,  Ω ≠  ℂ,
 of the complex plane can be mapped analytically, one-to-one and onto the unit disc. In 1851 as the theorem came to light in Riemann's innaugural dissertation address, [11], it stood at the threshold of tremendous new developments in complex analysis and geometry. Objections (by Weierstrass and others) to Riemann's proof  were  themselves to lead to whole new areas of mathematics, but a completely correct proof  of the mapping theorem would require another five decades.  
 
 That Riemann's statement of 1851 had been fundamentally correct was finally shown by Osgood, who in 1900 produced a proof which meets  our modern standards of rigor in analysis. See [8].  The proof of the theorem which has been standard in graduate courses and textbooks for about 70 years is, however, not Osgood's potential theoretic one but rather a function theoretic argument in which the mapping is shown, by the theory of normal families, to be the solution of an extremal problem. See [1] for example. In [13], Walsh credits Fejer and Riesz for the argument and Rado for its publication in 1923.  The standard proof is non-constructive because the mapping function is shown to be a subsequential limit who's existence is assured by the normal families argument. It is also out of reach for almost any undergraduate course in complex variables because of the time and background  required to introduce compactness in function space. See [10], [12] and [13] for more on the history of the theorem and its many connections with modern analysis.

A constructive approach

In the 1915 paper [7], Koebe gave a proof of the Riemann mapping theorem which constructed the conformal mapping by  the method of successive approximations introduced in [6]. In [9] Ostrowski gave another presentation of Koebe's proof  which with some minor variations appears in Bishop's "Foundations of Constructive Analysis", [2]. In the same paper, Ostrowski showed that the convergence of the algorithm
 in the proof would in general be asymptotically slow. The error in the  nth approximation to the mapping function would be on the order of 1/n. This estimate seems to have been taken as an indication of the 
 impracticality of the method and no numerical experiments were done for another fifty years.
 
 In 1979, E. Grassmann undertook an experimental study of Koebe's constructive proof in [3]. By using
 a method based on  three elementary mappings instead of only one, he observed the suprising phenomena 
 of fast initial convergence of the algorithm for a wide variety of domains. In [4] and [5], Henrici gives some mathematical foundation for these experimental conclusions. Grassmann's experiments were done with a Fortran program on the CDC 6400/6500/CYBER/174 at ETH in Zürich. One imagines that at the time, this was a powerful computing machine of a type not easily available to the general mathematical public.  The present software package running on a 64MB machine is quite up to the computing task and because Grassmann was meticulous in [3] in describing his algorithm, we are able to reproduce the results here. The combination of numerical and graphical capabilities makes for a nice classroom demonstration of  "what the Riemann mapping theorem says".  The proof of convergence of the method and therefore
of the Riemann mapping theorem can be achieved in an undergraduate course by an instructor who is
 inclined to plan ahead with that goal in mind. The required background results are Schwarz's lemma, Harnack's inequality, Hurwitz' theorem,  the general notion of a  continuous branch of the logarithm, and the Cauchy convergence criteria in the plane.

Contents

In the next section we will describe Grassmann's algorithm ( though not as carefully as he  did
 in [3]). In section 3 we  give the convergence proof assuming the above mentioned results. In section 4
 we give a Mathematica implementation of Grassmann's algorithm as described in [3]. (Bugs are of course
 my fault.) In section 5, some of Grassmann's examples and a few others are reproduced. We also give some warning to the potential experimenter about domains which will thwart the implementation.

2. Grassmann's Algorithm

Let 𝔻 denote the unit disc, \({z ∈ \(ℂ : \(\(|\)\(z\)\(|\)\(\(<\)\(1\)\)\)\)} . \ 
            We\ seek\ a\ conformal\ mapping\ 𝒻\ from\ an\ arbitrary\ simply\ connected\ domain\ \ Ω\ onto\ 𝔻 . \ It\ is\ enough\ to\ assume\ that\ Ω ⊂ 𝔻\),0 ∈ Ω,and that ∂Ω contains at least two points of ∂𝔻 .  An arbitrary domain can be transformed by a fractional linear transformation  T(z)=((az  +  b)/(cz  +  d)) to one satisfying these
\(assumptions . \ 
          In\ this\ section\ we\ will\ follow\ Grassmann\ in\ constructing\ a\ sequence\ of\ \ successively\ better\ approximations\ \ 𝒻\_n, \ 
      conformal\ on\ Ω, \ 
      which\ all\ satisfy\ \(𝒻\_n\) \((0)\) = 
        0, \ \(𝒻\_n\^′\) \((0)\) > 
          0\ and\ \ \(𝒻\_n\) \((Ω)\) ⊂ 𝔻 . \ 
            In\ section\ 3\ we\ will\ prove\ the\ convergence\ of\ the\ 𝒻\_n\ to\ the\ unique\ conformal\ mapping\ 𝒻\ from\ Ω\ onto\ 𝔻, \ 
      satisfying\ 𝒻 \((0)\) = 
        0\ and\ \(𝒻\^′\) \((0)\) > 
          0. \ The\ mapping\ \ 𝒻\_n\ is\ defined\ as\ a\ composition\ of\ elementary\ mappings\) 
\(𝒻\_n = \(𝒽\_\(\(n\)\(-\)\(1\)\(\ \)\)∘𝒽\_\(n - 2\)∘ ... \)∘𝒽\_0\ where\ each\ \ 𝒽\_j\ is\ conformal\ \ in\ a\ subdomain\ of\ 𝔻\ containing\ \(𝒻\_j\) \((Ω)\)\ and\ is\ normalized\ so\ that\) 
\(\(𝒽\_j\) \((0)\) = \(0\ and\ \(𝒽\_j\^′\) \((0)\) > 
          0. \ Each\ mapping\ \ \ 𝒽\_j\ is\ chosen\ to\ be\ one\ of\ three\ types\ depending\ on\ geometric\ criteria . 

Write\ Ω\_j\  = \(𝒻\_j\) \((Ω)\)\ \ and\ \ let\ \ z\_0\ denote\ a\ point\ of\ minimal\ modulus\ on\ ∂Ω\_j\  . \ 
            If\ there\ is\ a\ circle\ passing\ through\ \ z\_0\ which\ contains\)\)
\(Ω\_j\ , \ 
      find\ one\ of\ minimal\ radius . \ By\ our\ normalizations, \ 
      such\ a\ circle\ must\ intersect\ \ ∂𝔻\ and\ we\ choose\ 𝒽\_j\ to\ be\ the\ normalized\ conformal\ mapping\ from\ the\ intersection\ of\) 
\(the\ interior\ of\ this\ circle\ with\ 𝔻\ onto\ 𝔻 . \ 
        An\ example\ of\ this\ is\ shown\ in\ figure\ 1\ \ which\ is\ taken\ from\ the\ inverted\ ellipse\ computation\ in\ section\ 5. \ If\ no\ such\ circle\ can\ be\ found, \ we\ find\ a\ circle, 𝒞, \ 
    of\ maximal\ radius\ whose\ interior\ is\ disjoint\ from\ Ω\_j . \ If\ this\ circle\ intersects\ ∂𝔻\ we\ choose\ 𝒽\_j\ to\ be\ the\ normalized\ conformal\ mapping\ of\ the\ intersection\ of\ the\ \ exterior\ of\ 𝒞\ \ with\ 𝔻\ onto\ \ 𝔻 . \ 
        This\ is\ illustrated\ in\ figure\ 2\)
\(which\ is\ taken\ from\ the\ peanut\ example\ in\ section\ 5. \ \ In\ either\ of\ these\ two\ cases, \ 
    we\ say\ \ that\ the\ alogorithm\ has\ performed\ a\ circle\ \(\(routine\)\(.\)\(
\)\(
\)\)\)

Figure 1

Figure 2

If the circle 𝒞 of does not intersect ∂𝔻 we attempt to perform the ray routine which we will now describe. We check to see if the ray from (Subscript[z, 0]) through 
the center of 𝒞 intersects ∂(Subscript[Ω, j]). If not, then we choose  (Subscript[𝒽, j]) to be the conformal mapping from the intersection  of the complement of this ray with 𝔻 onto 𝔻. This is shown in figure 3 which is again taken from section 5.

Figure 3

\(In\ practice, \ 
    it\ is\ unwise\ to\ start\ the\ ray\ run\ at\ \ z\_0\ and\ in\ the\ implementation\ below\ we\ have\ let\ the\ ray\ start\ from\ \ the\ point\ \ P\ which\ is\ γ = \(\(3/
          4\)\(\ \)\(of\)\(\ \)\(the\)\(\ \)\(way\)\(\ \)\(from\)\(\ \)\(the\)\(\ \)\(center\)\(\ \)\(of\)\(\ \)\(𝒞\)\(\ \)\(to\)\(\ \ \)\(z\_0 . \ 
          Making\)\(\ \)\(γ\)\(\ \)\(closer\)\(\ \)\(to\)\(\ \)\(1\)\(\ \)\(increases\)\(\ \)\(the\)\(\ \)\(efficiency\)\(\ \)\(of\)\(\ \)\(the\)\(\ \)\(current\)\(\ \)\(ray\)\(\ \)\(routine\)\(\ \)\(but\)\(\ \)\(causes\)\(\ \)\(greater\)\(\ \)\(ditortion\)\(\ \)\(and\)\(\ \)\(decrease\)\(\ \)\(of\)\(\ \)\(smoothness\)\(\ \)\(of\)\(\ \)\(the\)\(\ \)\(image\)\(\ \)\(boundary . \ Decreasing\)\(\ \)\(γ\)\(\ \)\(decreases\)\(\ \)\(the\)\(\ \)\(efficiency\)\(\ \)\(of\)\(\ \)\(the\)\(\ \)\(current\)\(\ \ \)\(ray\)\(\ \)\(routine\)\(\ \)\(but\)\(\ \)\(may\)\(\ \)\(tend\)\(\ \)\(to\)\(\ \)\)\)
\(allow\ an\ earlier\ use\ of\ the\ circle\ routine\ in\ later\ iterations\ and\ \ thus\ increase\ the\ overall\ efficiency . \ 
        The\ line\ of\ the\ \ program\ in\ which\ this\ parameter\ γ\ is\ set\ to\ 3/4\ is\ in\ large\ blue\ \(\(type\)\(.\)\)\)

\(If\ the\ ray\ determined\ above\ intersects\ Ω\_j\), then we must resort to the function used by Koebe in his 1915 proof [7].  An example in which this is necessary is shown in 
\(figure\ 4. \ We\ refer\ to\ this\ case\ as\ the\ Koebe\ \(\(routine\)\(.\)\)\)

Figure 4

The Koebe Routine

\(In\ the\ Koebe\ routine, \ 
      the\ point\ z\_0\ is\ moved\ to\ the\ origin\ by\ a\ fractional\ linear\ transformation\ \ S\ which\ preserves\ the\ unit\ disc . \ 
            The\ domain\ Ω\_j\ is\ transformed\ to\ a\ simply\ connected\ domain\ on\ which\ we\ can\ define\ a\ continuous\ branch\ of\ the\ square\ root . \ 
            Applying\ this\ square\ root\ mapping\ opens\ up\ the\ domain\ near\ the\ origin\ and\ then\ \ another\ disc\ preserving\ fractional\ linear\ transformation\ is\ applied\ to\ force\ the\ normalizations\ \(h\_j\) \((0)\) = 0\ and\ \(h\_j\^′\) \((0)\) > 
          0. \ In\ the\ discretized\ algorithm\ it\ is\ actually\ the\ point\ P\ that\ we\ deal\ with\ rather\ than\ z\_0 . \ 
              The\ algorithm\ first\ rotates\ the\ whole\ picture\ so\ that\ P\ can\ assumed\ to\ be\ real\ and\ negative . \ The\ mapping\ S\ is\ then\) 
                                                                                          S (z)=((z - P)/(1 - Pz)),
and the  Koebe routine mapping is 
                                                                           \(\(h\_j\) \((z)\) = \(\@\(S \((z)\)\) - \@\(\(|\)\(P\)\(|\)\)\)\/\(1\  - \(\@\(\(|\)\(P\)\(|\)\)\) \@\(S \((z)\)\)\)\) ,
 where the square root is our continuously defined branch .  The whole picture is then rotated back to the correct position .   The definition of a continuous branch of the argument is shown in green type in the program listing.

A  Circle routine map

\(To\ write\ down\ the\ mapping\ from\ the\ circle\ routine, \ 
    send\ one\ intersection\ point\ to\ the\ origin\ and\ the\ other\ to\ the\ point\ at\ infinity\ by\ a\ fractional\ linear\ transformation . \ 
        By\ a\ rotation\ and\ power\ function, \ 
    map\ the\ resulting\ angular\ wedge\ to\ the\ upper\ halp\ plane . \ 
        Map\ to\ the\ unit\ disk\ \ and\ fix\ the\ normalizations\ by\ a\ further\ fractional\ linear\ transformation . \ 
        In\ the\ numerical\ routine, \ 
    the\ picture\ is\ rotated\ to\ the\ convenient\ position\ shown\ below . \ Denote\ the\ angle\ interior\ to\ the\ domain\ between\ the\ two\ circles\ by\ ϕ, \ 
    and\ the\ upper\ intersection\ by\ Z . \ 
        Then\ the\ desired\ mapping\ \(h\_j\) 
      is\ given\ by\)
                                                                      \(T \((z)\) = \ \(\(z - Z\)\/\(z - Z\&_\)\) . \(\(z\_0 - 
                    Z\&_\)\/\(z\_0 - Z\)\)\),
                                               (Subscript[h, j]) (z)= \(\(T \((z)\)\^\(π\/ϕ\) - 
                      T \((0)\)\^\(π\/ϕ\)\)\/\(T \((z)\)\^\(π\/ϕ\) - \(T \((0)\)\^\(π\/ϕ\)\)\&_\)\). \(\(T \((1)\)\^\(π\/ϕ\) - \(T \((0)\)\^\(π\/ϕ\)\)\&_\)\/\(T \((1)\)\^\(π\/ϕ\) - 
                      T \((0)\)\^\(π\/ϕ\)\)\).
The figure depicts curves corresponding to circles and rays centered at the origin in the unit disc.

generating the figure

The figure takes some time for the computer to draw. So the cells are switched off.

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]

<< "ProgrammingInMathematica`ComplexMap`"

x =  ((- .9  +   .3) (- .9  -   .3)  + 
            1)/(2 (- .9))
y =  Sqrt[1  -  x^2]
theta  = N[-ArcSin[y]  +  ArcSin [y/ .3], 16]
a = (- .6  -  (x -  I y))/(- .6  -  (x  +  
                I y))
Z1 =  x  +  I y
Tinverse[
        w_] :=  ((w/a) Conjugate[Z1]  -  Z1)/((w/a) - 
            1)
T[0] = a Z1/Conjugate[Z1]
T[1] =  a (1  -  Z1)/Conjugate[1 - Z1]
b =  (T[1]^(Pi/theta)  -  
            Conjugate[T[0]]^(Pi/theta))/(T[1]^(Pi/theta)  -  
            T[0]^(Pi/theta))
Pinverse[
        w_] :=  (((w/b) Conjugate[T[0]]^(Pi/theta)  - 
                T[0]^(Pi/theta))/((w/b)  - 1))^(theta/
            Pi)
inverseCrescent[z_] := Tinverse[Pinverse[z]]

graph0 = PolarPlot[1, {t, 0, 2 N[π]}];
graph1 = Show[Graphics[Circle[{- .9, 0},  .3]]];
graph2 = 
        PolarMap[(inverseCrescent[#]) &, {0,  .9,  .1}, {-Pi,  Pi, 
            2 Pi/30}];
graph3 = 
        PolarMap[(inverseCrescent[#]) &, { .9,  .99,  .01}, {-Pi,  
            Pi, 2 Pi/30}];
graph5 = 
        PolarMap[(inverseCrescent[#]) &, { .999,  .9999,  .0001}, {-Pi,  Pi, 2 Pi/30}];

Show[graph0,graph1,graph2,graph3,PlotRange->{-1,1}]

A ray routine Map

\(For\ the\ ray\ routine, \ 
    we\ again\ rotate\ the\ picture\ for\ convenience\ as\ shown\ below . \ 
        Let\ 𝓉\ denote\ the\ tip\ of\ the\ slit\ in\ the\ interior\ of\ the\ unit\ circle . \ 
        The\ function\ \((usually\ named\ after\ Koebe\ but\ not\ to\ be\ confused\ with\ the\ function\ from\ the\ Koebe\ routine)\)\)
                                                                               k (z) =  (z/(((1 - z))^2))
\(maps\ the\ unit\ disc\ conformally\ onto\ the\ plane\ minus\ the\ real\ interval\ [\(-∞\), \(-\(1\/4\)\)] . \ The\ inverse\ function\ is\)
                                                                   \(\(k\^\(-1\)\) \((z)\) = \(\(2  z\  + \ 1\  - \ \@\(4  z\  + \ 1\)\)\/\(2  
                  z\)\) . 
The\ required\ mapping\ \ is\ then\)
                                                                       (Subscript[h, j]) (z)=  (k^(-1))(((((t + 1))^2)/(4  t)) k (z)).
Curves  corresponding to circles and rays centerd at the origin in the unit disc are shown in the figure below.

Generating the figure.

Again, the figure may take time to be redrawn. The cells are switched off.

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]

<< "ProgrammingInMathematica`ComplexMap`"

shiftedarg[z_] :=  
        If [(-Pi/2) <  Arg[z]  <= Pi,  Arg[z],  
          Arg[z]  +  2  Pi];
invpower[z_] :=  
        Exp[(theta/Pi) (Log[Abs[z]]  +  
                I  shiftedarg[z])];
k[z_] := (z/(((1 - z))^2));
kinverse[z_] := ((2 z + 1 - (Sqrt[4 z + 1]))/(2 z));
S[z_, t_] := ((z - t)/(1 - t z));
smap[z_] := 
        kinverse[(((((Abs[- .6] + 1))^2) k[z])/(4 Abs[- .6]))];
sinverse[w_] :=  
        kinverse[(4 Abs[- .6] k[
                  w])/(((Abs[- .6] + 1))^2)];

graph0 = PolarPlot[1, {t, 0, 2 N[π]}]

graph1 = 
      PolarMap[(sinverse[#]) &, {0.1, 1,  .08}, {-Pi  ,  Pi, 
          2 N[Pi]/60}]

Show[graph0,graph1,Graphics[{AbsoluteThickness[2],Line[{{-1,0},{-.6,0}}]}]]

3. The Convergence Proof

The proof.

We follow the presentations in [3] and [5].

The convergence of the discrete version of the algorithm is clear. It requires at most a number of iterations equal to the number of discretized boundary points. 

To prove that the continuous version of the algorithm produces a convergent sequence we keep in mind the 
notation and normalizations from the previous section.
Let (Subscript[Ω, n]) denote the result of n iterations applied to the original domain Ω. Let (Subscript[a, n]) denote the minimum modulus of points on ∂(Subscript[Ω, n]). Let (Subscript[h, n]) denote the transformation to be performed next and let  (Subscript[f, n])= (Subscript[h, n - 1])∘...∘(Subscript[h, 0]). For each of the three types of mappings, the inverse functions (Subscript[h, n]^(-1)) extend analytically to the whole unit disc. The extension is not necessarily one-to -one, but maps the unit disc into itself. Applying Schwarz's lemma to the function (Subscript[h, n]^(-1)) we find that 
|(Subscript[h, n])(z)| > |z| and that  |((Subscript[h, n])^′)(0)|>1. By the normalization on the mappings in the algorithm we actually have ((Subscript[h, n])^′)(0)>1. The sequences {(Subscript[a, n])} and {((Subscript[f, n])^′)(0)} are thus increasing and bounded respectively by 1 and by (1/(Subscript[a, 0])) (by Schwarz's lemma applied to (Subscript[f, n])((Subscript[a, 0]) z)).  Denote the limits respectively by 𝒶 and φ. Though the function (Subscript[h, n]) need not be defined at a given  point of ∂(Subscript[Ω, n]), the analyticity of  (Subscript[h, n]^(-1)) in 𝔻 implies that  {(Subscript[a, n])} is in fact strictly increasing.

We claim that 𝒶=1. If 𝒶< 1, consider the set of simply connected sub-regions Ω of the unit disc containing the origin with 

                                                                     (a/2)≤min{|z|:z ϵ ∂Ω}≤𝒶.  
 
 Let Φ denote the set of all transformations which
might be used by the algorithm's first step on any such sub-region. By Schwarz's lemma, 

                                                                      σ= inf{(h^′)(0):h ϵ Φ}>1.

 But for a  sufficiently large n we have (Subscript[a, n])≥𝒶/2 and  ((Subscript[f, n])^′)(0)> φ/σ which implies ((Subscript[f, n + 1])^′)(0)=((Subscript[h, n])^′)(0)((Subscript[f, n])^′)(0)>φ. 
 This is a contradiction and shows that we must have 𝒶=1.

We have shown that the domains  (Subscript[Ω, n]) tend to the unit disc and now will show that on a given compact subset K⊂Ω the functions (Subscript[f, n]) converge uniformly to 𝒻. If n>m then the function ((Subscript[f, n])/(Subscript[f, m])) is analytic and not zero in Ω since both the numerator and denominator
have exactly one simple zero at the origin. By Schwarz's lemma we have |((Subscript[f, n])/(Subscript[f, m]))|>1 in Ω. Because Ω is simply connected there is a unique analytic branch of 

                                                                           p(z) =log(((Subscript[f, n])(z))/((Subscript[f, m])(z))) 
                                                           
 defined in Ω which agrees with the principle value 
 
                                                                             p(0) =Log ((((Subscript[f, n])^′)(
            0))/(((Subscript[f, m])^′)(0))) >0 
  
 at the origin and we have
 
                                                                          Re p(z) = log|(((Subscript[f, n])(z))/((Subscript[f, m])(z)))| >0 
  
  in Ω.
The function (Subscript[f, n])∘(Subscript[f, m]^(-1)) maps a disk of radius (Subscript[a, m+1 ])onto the unit disc while fixing the origin so another 
application of Schwarz's lemma shows that 

                                                                               \(TraditionalForm\`\(\(f\_n\^′\)(0)\)\/\(\(\ \ \ \ \ \)\(\(f\_m\^′\)(0)\)\)\)< (1/(Subscript[a, m + 1]))<(1/(Subscript[a, m])). 

It follows that p(0) < Log(1/(Subscript[a, m])) and therefore by Harnack's inequality for positive harmonic functions that there is a positive constant c=c(K) such that

                                                                                   |p(z)| ≤cLog(1/(Subscript[a, m]))
 
 for all z in K.

The elementary inequality

                                                                                     |(e^u)-1|≤(e^(|u|))|u| 
 
 follows from an integration and shows that
 
                                                                                |(((Subscript[f, n])(z))/((Subscript[f, m])(z)))-1|≤\(TraditionalForm\`c\/a\_m\%c\)Log(1/(Subscript[a, m]))  
  
  for all z in K. Another elementary inequality,
  
                                                                                         Log(1/x)≤(1/x)-1 
   
   valid for all x>0, together with |(Subscript[f, m])(z)|<1  then shows that
   
                                                                              |(Subscript[f, n])(z)-(Subscript[f, m])(z)| ≤\(TraditionalForm\`c\/a\_m\%\(c + 1\)\)(1-(Subscript[a, m])) 
    
    for z in K.

Since (Subscript[a, m])→1, the Cauchy criteria gives the uniform convergence of the (Subscript[f, n]) on K to a limit function f
which is one-to one by an application of Hurwitz's theorem, maps Ω onto the unit disc, and has the correct normalization. The proof is complete.

Some thoughts on teaching the proof.

Most one semester undergraduate courses in complex analysis probably cover the maximum principle
for  analytic and harmonic functions, so the inclusion of Schwarz's lemma is not a problem.  By 
the end of the semester, the students  may have also  grasped the concept of a general continuous branch of the 
logarithm and the general notion of simple connectivity can be taken on faith. The elementary mappings used in Grassmann's algorithm can be covered as examples during the appropriate standard part of the course although the Koebe routine presents the difficulty of a generalized branch of the square root. The tough part is to include both Harnack's principal, requiring the Poisson integral formula  for harmonic functions and Hurwitz's theorem which is  a consequence of the argument principal. Assuming one plans 
carefully to cover all these results in time to present the Riemann Mapping theorem, will students 
get the punch line? Hopefully the program and demonstration examples in the next two sections
may be of some help. But even without all the background, a demonstration of Grassmann's experiments
may be useful to give a general impression of what the Riemann Mapping theorem says.

4. The Program

The following program is freely available at 
http: // math . claremontmckenna . 
                      edu/moneill/.  It should run  on smooth boundaries 
\(which\ are\ not\ exactly\ circular\ or\ piecewise\ circular . \ Make\ sure\ that\ the\ data\ entered\ as\)
\(pointlist\ \*"""" 
        is\ a\ list\ of\ cartesian\ coordinates\ {x, y} . \ 
          The\ points\ must\ lie\ in\ the\ unit\ disk\ and\ at\ least\ two\ of\ them\ must\ have\)
\(\(\ \)\(modulus\ 1. \ Also\ make\ sure\ that\ there\ are\ no\ repeated\ points\ anywhere\ in\ the\ list . \ \ Non - 
        smooth\ or\ very\ complicated\)\)
 boundaries  may cause  problems with the functions  circleq, \(slitq\ or\ with\ the\ Koebe\ Routine\ through\ the\ function\) 
Newarg or the built in Arg function .     Most problems are related to  trying to evaluate a branch of
\(\(\ \)\(the\ argument, \ 
      on\ a\ discretized\ boundary\  . \ 
\ 
          The\ package\ is\ meant\ for\ demonstration\ purposes\ only . \ 
          It\ might\ be\ fun\ to\ construct\ a\ version\ which\ works\)\(\ \)\)
on doubly connected domains. See  [5].

BeginPackage["Osculation`"]

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]


GRASSMANN::usage = "GRASSMANN[pointlist,tol,maxiterates] \
produces a graphical display
		of successive approximations, obtained by the Grassman version of Koebe's \
osculation algorithm, to the Riemann mapping of a simply connected domain \
onto
		the unit disk. The boundary of the domain  is represented as a list of \
cartesian coordinates which are input as
		the variable pointlist. There should be no repetition in the list of \
coordinates. It is asumed that the domain is contained in the unit disk and \
that  at least two points in the 
			discretization of
			 the boundary lie on the unit circle. The algorithm halts if the  number \
of iterations exceeds the variable maxiterates or if the
			minimal modulus of the discretized boundary exceeds the variable tol which \
must
			be a real number between 0 and 1.";


Begin["`Private`"]
   
				 
		GRASSMANN[pointlist_, tol_, iteratenumber_]:=Module[{g1, unitcircle, points, complexpoints, iterations, 
              minmod, appRiemann, S, multiplicationfactor, function0, 
              circleq, slitq, newArg, koebesqrt, k, kinverse, modulisquared, 
              minpairs, zi0, zi1, pointsminusaminpair, Ajlist, A, r, c, 
              thirdpoint, testpoint, mfactor, P, x, y, Z1, T, crescent, 
              crescentappliedtorotatedpoints, Kroutineappliedtorotatedpoints, 
              smap, smapappliedtorotatedpoints},
		
		  
		
		points=pointlist;
		
	  
			
		unitcircle = 
                PolarPlot[1, {t, 0, 2 N[π]}, 
                  DisplayFunction -> Identity];
		
		\(g1 = 
                Show[Graphics[Map[Point, Append[points, First[points]]]], 
                  unitcircle, 
                  AspectRatio -> 
                    Automatic, 
				DisplayFunction\  -> $DisplayFunction]\);
		
		  
		
		complexpoints = (#1[[1]] + I #1[[2]] &) /@ 
                  points;
		
		  
		
		iterations = 0;
		
		minmod = 0;
		
		
		  
		
		S[z_, t_]:=((z - t)/(1 - t z));
		
		
		  
		
		multiplicationfactor[pt_, cntr_, radius_] := 
                If[Abs[pt - cntr] >= 
                    radius, -Conjugate[cntr]/Abs[cntr], 
                  Conjugate[cntr]/Abs[cntr]];
		
	  
		
		\(circleq[h_, k_, s_] := 
                h\^2 + k\^2 > 1 || 
                  0 < h\^2 + k\^2 < 1 && s > 1 - \((h\^2 + k\^2)\)\^ .5\);
		
		  
			
				slitq[complexlist_, a_] := 
                Module[{args, jumps, 
                    badjumppositions}, 

args = 
                    Arg /@ Append[complexlist, 
                        First[complexlist]]; 
jumps = (Abs[#1[[1]] - #1[[2]]] &) /@ 
                      Partition[args, 2, 
                        1]; 
badjumppositions = 
                    complexlist[[
                      Flatten[
                        Position[(#1 >  .5 &) /@ jumps, 
                          True]]]]; 

Min[Abs /@ badjumppositions] > 
                    Abs[a]

                         ];
					  
		
		newArg[ls_] := 
                  Module[{contarg}, 
                    contarg[k_] := 
                      If[k == 1, 
                        Arg[ls[[1]]], 
                        contarg[k] = 
                          If[Abs[
                                Arg[ls[[k]]] - Arg[ls[[k - 
                                        1]]]] <  .5, 
                            contarg[k - 1] + 
                              Arg[ls[[k]]] - Arg[ls[[k - 1]]], 
                            If[Arg[
                                    ls[[k]]] - Arg[ls[[k - 
                                        1]]] ≥  .5, contarg[k - 1] + 
                                Arg[ls[[k]]] - Arg[ls[[k - 
                                      1]]] - 2 N[π], 
                              contarg[k - 1] + 
                                Arg[ls[[
                                    k]]] - 
                                Arg[ls[[
                                    k - 1]]] + 
                                2 N[π]]]]; 

Table[
                      contarg[n], {n, 1, 
                        Length[
                          ls]}]

                          ];
		
		  
		koebesqrt[ls_, 
                  n_] := ((Sqrt[Abs])[#1] &) /@ ls (Exp[((I #1)/2)] &) /@ 
                    RotateRight[newArg[RotateLeft[ls, n]], n];
		
		k[z_] := (z/(((1 - z))^2));
		kinverse[
                  z_] := ((2 z + 1 - (Sqrt[4 z + 1]))/(2 z));
		  
		
		While[iterations < iteratenumber && minmod < tol,
			
			  
				
			modulisquared = (#1[[1]] #1[[1]] + #1[[2]] #1[[2]] &) /@ 
                        Join[{points[[-1]]}, points];
			
			  
				
			minpairs=\(Transpose[{\(Join[{points[[\(-1\)]]}, points]\)[[Flatten[
                              Position[modulisquared, 
                                Min[modulisquared]]]]], \(Join[{points[[\(-1\)]]}, 
                              points]\)[[
							If[
								MatchQ[
									Flatten[
                                  Position[modulisquared, 
                                    Min[modulisquared]]], {1, ___}
								], 
								
									Flatten[
                                Prepend[
										Rest[
											Position[modulisquared, Min[
                                        modulisquared]]
										] - 1, Length[modulisquared] - 
                                    1]], 
								
								Flatten[Position[modulisquared, Min[modulisquared]]] - 
                                1
							]
								
								]]}]\);
			
			
			  
			
			
			zi0 = (#1[[1]] + I #1[[2]] &)[
                        minpairs[[1, 
                          1]]];
			

			
			zi1 = (#1[[1]] + I #1[[2]] &)[
                        minpairs[[1, 
                          2]]];
			
		  
			
			pointsminusaminpair=Delete[
                        points, {Flatten[
                            Position[points, 
                              minpairs[[1, 
                                1]]]], 
                          Flatten[
                            Position[points, 
                              minpairs[[1, 
                                2]]]]}];
		
			complexpointsminusaminpair  =  
                      Map[(#1[[1]] + 
                              I #1[[2]] &), pointsminusaminpair];	
			
				functionA[
                        z_] := ((z - zi0)/(z - 
                                zi1)) (Abs[(z - zi1)/(z - 
                                  zi0)]);
				
						
			Ajlist =  
                      Map[functionA, complexpointsminusaminpair];
			
			
			
			A = Ajlist[[
                          Flatten[
                            Position[Im /@ Ajlist, 
                              Min[Im /@ 
                                  Ajlist]]]]][[1]];
			
			

			
			
			r = (( .5 Abs[zi0 - zi1])/Abs)[Im[A]];             
			
		
			c =  .5 (zi0 + 
                            zi1 + ((I (zi0 - zi1) Re[A])/Im)[A]);         
			
	
			   
		
			testpoint = Complement[{pointsminusaminpair[[1]], 
                            pointsminusaminpair[[2]]}, pointsminusaminpair[[
                            Flatten[
                              Position[Im /@ Ajlist, 
                                Min[Im /@ 
                                    Ajlist]]]]]][[1]];
			
		
			
			mfactor = 
                      multiplicationfactor[(#1[[1]] + I #1[[2]] &)[
                          testpoint], c, r];
		
			
			
			rotatedcomplexpoints = 
                      multiplicationfactor[(#1[[1]] + I #1[[2]] &)[
                            testpoint], c, r] complexpoints;
			
			  
			
			P = mfactor c +  .75 r;
			
				  									
			
			If[circleq[Re[c], Im[c], r],
				
				   
				
				If[r >  10^6,
					
					lineplot  =  
                                ImplicitPlot[(Re[zi0]  -  
                                        Re[zi1]) q  -  (Im[ 
                                        zi0]  -  Im[ zi1]) p  -  
                                      Im[ zi0] (Re[ zi0]  -  
                                        Re[ zi1])  +  
                                      Re[ zi0] (Im[ zi0]  -  
                                        Im[zi1]) == 0, {p, -1.2, 
                                    1.2}, {q, -1.2, 1.2},  
                                  DisplayFunction  -> Identity];
					
						x = Re[mfactor  zi0];
						
						y = Sqrt[1 - x^2];
						
						
						Z1 = x  +  I  y;
						
					\(phi = 	ArcCos[x]\);
						
						\(Show[lineplot, g1, 
                                unitcircle, AspectRatio -> Automatic, 
                                Axes -> True, 
                                PlotRange -> {{\(-1.2\), 
                                      1.2}, {\(-1.2\), 1.2}}, 
                                DisplayFunction
					 -> $DisplayFunction]\);
						
					  	
					
				T[
                                  z_] := (((z - Z1) (mfactor zi0 - Conjugate[Z1]))/((z - Conjugate[Z1]) (mfactor zi0 - Z1)));
				
				crescent[
                                  z_] := ((T[z]^(N[π]/phi) - 
                                        T[0]^(N[π]/phi)) (T[1]^(N[π]/phi) - 
                                        Conjugate[
                                        T[0]^(N[π]/phi)]))/((T[z]^(N[π]/phi) - 
                                        Conjugate[
                                        T[0]^(N[π]/phi)]) (T[1]^(N[π]/phi) - T[0]^(N[π]/phi)));
				
				crescentappliedtorotatedpoints = 
                                crescent /@ rotatedcomplexpoints;
				
					  
				complexpoints = 
                                (crescentappliedtorotatedpoints/mfactor);
				
				points = ({Re[#1], 
                                        Im[#1]} &) /@ complexpoints;
				
				iterations = iterations + 1;
				
				minmod = Abs[zi0];
				
					  
																																																				
		\(g1 = 
                                Show[Graphics[
                                    Map[Point, 
                                      Append[\(({Re[#1], 
                                        Im[#1]} &)\) /@ \(crescentappliedtorotatedpoints\/mfactor\), 

                                        First[\(({Re[#1], 
                                        Im[#1]} &)\) /@ \(crescentappliedtorotatedpoints/mfactor\)]]
								]], unitcircle, 
                                  AspectRatio -> 
                                    Automatic, 
						PlotRange\  -> {{\(-1\), 1}, {\(-1\), 1}}, DisplayFunction\  -> $DisplayFunction]\),
						
			  			
						
						 \(Show[
                                  Graphics[Circle[{Re[c], Im[c]}, r]], g1, 
                                  unitcircle, AspectRatio -> Automatic, 
                                  Axes -> True, 
                                  PlotRange -> {{\(-1\), 1}, {\(-1\), 
                                        1}}, 
                                  DisplayFunction
					 -> $DisplayFunction]\);
				
				x = 
                                Re[(((mfactor  c + r) (mfactor  c - r) + 1)/(2  mfactor  c))];
			
				y = (Sqrt[1 - (x^2)]);
				
				Z1 = x + I y;
				
			
				\(phi = 
                                If[Positive[Re[mfactor\ c]], 
                                  N[π] - \((ArcSin[y] - ArcSin[y\/r])\), 
                                  If[
							Re[
                                        mfactor\ c]\  <= \(-1\), 
                                    N[π] - \((ArcSin[y] + 
                                        ArcSin[y\/r])\), 
                                    If[Re[
                                        mfactor\ c\ ] <= \(-Sqrt[1 - r^2]\), 
                                      N[π] - \((ArcSin[y] + 
                                        ArcSin[y\/r])\), \ 
                                      ArcSin[y/r]\  - \ ArcSin[y]]]]\);
		
				T[
                                  z_] := (((z - Z1) (mfactor zi0 - Conjugate[Z1]))/((z - Conjugate[Z1]) (mfactor zi0 - Z1)));
				a = T[0]^(N[π]/phi);
					b = T[1]^(N[π]/phi);
					c = Conjugate[a];
				pow[z_] :=  
                                Exp[N[π]/phi Log[((z - 
                                        Z1) (mfactor zi0 - 
                                        Conjugate[Z1]))/((z - 
                                        Conjugate[Z1]) (mfactor zi0 - 
                                        Z1))]];	
				crescent[
                                  z_] := ((pow[z] - a) (b - 
                                        c))/((pow[z] - c) (b - 
                                        a));
				
				crescentappliedtorotatedpoints = crescent /@ rotatedcomplexpoints;
				
				complexpoints = 
                                (crescentappliedtorotatedpoints/mfactor);
				
				points = ({Re[#1], 
                                        Im[#1]} &) /@ complexpoints;
				
				iterations = iterations + 1;
				
				minmod = Abs[zi0];
										
																																												
				\(g1 = 
                                Show[Graphics[
                                    Map[Point, 
                                      Append[\(({Re[#1], 
                                        Im[#1]} &)\) /@ \(crescentappliedtorotatedpoints\/mfactor\), 

                                        First[\(({Re[#1], 
                                        Im[#1]} &)\) /@ \(crescentappliedtorotatedpoints/mfactor\)]]
								]], unitcircle, 
                                  AspectRatio -> 
                                    Automatic, 
						PlotRange\  -> {{\(-1\), 1}, {\(-1\), 1}}, DisplayFunction\  -> $DisplayFunction]\)],
				
				
				
		  
				
				
				
				If[slitq[rotatedcomplexpoints, 
                              mfactor c +  .75 r],Kroutineappliedtorotatedpoints = (S[#1, (Sqrt[Abs])[P]] &) /@ koebesqrt[(S[#1, P] &) /@ 
                                      rotatedcomplexpoints, Flatten[
                                        Position[modulisquared, 
                                        Min[modulisquared]]][[1]]];
					
					  
					
					complexpoints=(Kroutineappliedtorotatedpoints/mfactor);
					
					points = ({Re[#1], 
                                        Im[#1]} &) /@ complexpoints;
					
					iterations = iterations + 1;
					
					minmod = Abs[zi0];
					
					
				  
					
					
					g1 = 
                                Show[Graphics[
                                    Map[Point, 
                                      Append[points, First[points]]]], 
                                  unitcircle, AspectRatio -> Automatic, 
                                  DisplayFunction -> $DisplayFunction],
					
					
					
						  
					
					Show[g1, 
                                Graphics[
                                  Line[{{Cos[Arg[c]], 
                                        Sin[Arg[c]]}, {Abs[P] Cos[Arg[c]], 
                                        Abs[P] Sin[Arg[c]]}}]]];
					
					
					
				
					
					smap[z_] := 
                                kinverse[(((((Abs[P] + 1))^2) k[z])/(4 Abs[P]))];
					
					smapappliedtorotatedpoints = 
                                smap /@ rotatedcomplexpoints;
					
					   
					
					complexpoints = 
                                (smapappliedtorotatedpoints/mfactor);
					
					points = ({Re[#1], 
                                        Im[#1]} &) /@ complexpoints;
					
					iterations = iterations + 1;
					
					minmod = Abs[zi0];
					
					  
					
					g1 = 
                                Show[Graphics[
                                    Map[Point, 
                                      Append[points, First[points]]]], 
                                  unitcircle, AspectRatio -> Automatic, 
                                  DisplayFunction -> $DisplayFunction]]]];
		
		
		   
		
		
		minmod]
End[]
EndPackage[]

5. Examples

Though we have proven that the algorithm works perfectly well in its continuous version, discretization introduces significant difficulties. This can be seen in the last three examples which are, nevertheless, useful for giving a good idea of how Koebe's proof works. To see the examples, just evaluate the cells in order.

A Peanut

For this class of examples, the use of the circle routine is immediate and therefore the fast initial convergence is evident.

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]

g1 = 
        PolarPlot[1 -  .5  Cos[t]^2, {t, 0, 2 N[π]}, 
          PlotPoints -> 1000,  
          DisplayFunction -> Identity];

points = 
        Drop[FullForm[g1][[1, 1, 1, 1, 
            1]], -1];

ListPlot[points, AspectRatio -> Automatic]

GRASSMANN[points,  .95, 10]

A Polygon

The circle routine is not thrown off by a straight piece of boundary. (At least not in this example.)

\(\(pointlist = 
        Join[Drop[
            Table[\((1 - t)\) {Cos[N[Pi]/12], \ Sin[N[Pi]/12]}\  + \ 
                t {Cos[N[Pi]/1.8], Sin[N[Pi]/1.8]}, {t, 0, 1, 
                1/500}], \(-1\)]\ , 
		Drop[
            Table[\((1 - t)\)\ {Cos[N[Pi]/1.8], 
                    Sin[N[Pi]/1.8]}\  + 
					\ t\ {\(-Cos[
                        N[Pi]/12]\), \ Sin[N[Pi]/12]}, {t, 0, 1, 
                1/500}], \(-1\)], 
		Drop[
            Table[\((1 - t)\) {\(-Cos[N[Pi]/12]\), \ Sin[N[Pi]/12]}\  + \ 
                t {\(-Cos[N[Pi]/12]\), \(-Sin[N[Pi]/12]\)}, {t, 0, 1, 
                1/500}], \(-1\)]\ , 

          Drop[Table[\((1 - 
                      t)\) {\(-Cos[N[Pi]/12]\), \ \(-Sin[N[Pi]/12]\)}\  + \ 
                t {Cos[N[Pi]/12], \(-Sin[N[Pi]/12]\)}, {t, 0, 1, 
                1/500}], \(-1\)], \ 

          Drop[Table[\((1 - t)\) {Cos[N[Pi]/12], \ \(-Sin[N[Pi]/12]\)}\  + \ 
                t {Cos[N[Pi]/12], Sin[N[Pi]/12]}, {t, 0, 1, 
                1/500}], \(-1\)]];\)\)

ListPlot[pointlist, AspectRatio -> Automatic]

GRASSMANN[pointlist,  .95, 10]

An ellipse

Although the algorithm appears to make short work of the elongated ellipse, computation of the inverse mapping  would pose great numerical difficulties because of the size of the derivates of the circle maps near the ends of the ellipse.

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]

g1 = 
        PolarPlot[(1  +  24 ( Sin[t])^2)^(- .5), {t, 0, 
            2 N[π]}, PlotPoints -> 1000];

points = 
        Drop[FullForm[g1][[1, 1, 1, 1, 
            1]], -1];

GRASSMANN[points,  .95, 5]

An inverted ellipse

A five to one ellipse is inverted in the unit circle.
 Again, the algorithm converges rapidly at first.
  The computaton of the inverse would pose problems here near the sharp inverted tips.

Needs["Graphics`Graphics`"]
Needs["Graphics`ImplicitPlot`"]

g1 = 
        PolarPlot[ .2 (1  +  24 ( Sin[t])^2)^( .5), {t, 0, 
            2 N[π]}, PlotPoints -> 1000];

points = 
        Drop[FullForm[g1][[1, 1, 1, 1, 
            1]], -1];

GRASSMANN[points,  .95, 8]

A Daisy

This asymmetric daisy, though it is a starlike domain, forces the algorithm to use the Koebe routine several times before being able to switch to the ray routine and then finally the circle routine. Grassmann remarks that the example was chosen in "an only partially successful attempt to make the algorithm fail in practice". For a better picture but more computing time change the increment
(.0001) in pointlist1 to (.00001). With the current increment, the huge distortion near the inward pointy tips of the daisy are apparent. Grassmann notes in his paper that the use of the Koebe routine makes the computation of the inverse to the approximation numerically problematic. The customized choice of the branch of the logarithm together with a small round off error can lead to gigantic errors in the computation of the inverse. 

Changing 5 to 7  (in both positions) causes the program to fail quickly. And  the more so if the number of points is increased. The performance is better with 10 petals and increment .00001 for t. I believe that is because there is some luck involved in the choice of the initial boundary point of minimal modulus. The problems are arising for at least two reasons. One reason is  tha tthe very sharp inward pointing cusps are causing our circle 𝒞 from the circle routine to have zero radius (or at least practically zero). In some cases this is  combining with round off error in the branch of the argument to defeat the slitq criteria and fool the routine into doing a slit mapping when it should do the Koebe map. Another cause of the catastrophic errors is  that the very sharp cusps are making rounding errors dominant in computing the value of the argument in the branch of the squareroot for the Koebe map.

pointlist1 = 
        Drop[Table[{( .2 +  .8 Abs[Sin[5 N[π] t]]) Cos[
                  2 N[π] (t^1.2)], ( .2 +  .8 Abs[
                        Sin[5 N[π] t]]) Sin[
                  2 N[π] (t^1.2)]}, {t, 0, 
              1,  .0001}], -1];

ListPlot[pointlist1, AspectRatio -> Automatic]

GRASSMANN[pointlist1,.95,8]

Another problem example.

A domain bounded by circular arcs gives the case testing in the program a hard time.
(Catastrophically hard.) Similarly, the  discretized ray routine would have a difficult (or impossible) time mapping a slit disc.

\(z[s_] := \ 
      Which[0 <= Mod[s, \ 2\ Pi] < Pi, Cos[s]\  + \ I\ Sin[s], \ 
          Pi <= Mod[\ s\ , 2\ Pi] < \ 
            5\ \ Pi\ /4\ , \ \(-\((3/4)\)\)\  - \ \((1/4)\)\ \((\ 
                Cos[4\ s\ ]\  + \ I\ Sin[\ 4\ s])\), 
		\((5\ Pi\ /
                  4)\)\  <= \ Mod[s, 2\ Pi]\  < \ 
            7\ Pi\ /4, \ \((I/2)\) \((\ \(-\ Cos[2\ s\ ]\)\  + \ 
                I\ Sin\ [\ 2\ s])\), \ \((7\ Pi\ /4)\)\  <= 
            Mod[\ s, \ 2\ Pi]\  < \ 
            2\ Pi, \ \((3/
                4)\)\  + \ \((1/4)\)\ \((Cos[4\ s]\  + \ 
                  I\ Sin[4\ s])\)] - \  .75  I\)

angle =  Table[2  N[Pi] k/1000, {k, 1, 1000}];

pointlist = 
        Map[({Re[#], Im[#]}) &, 
          Map[z, angle]/(Max[Abs[Map[z, angle]]])];

ListPlot[pointlist]

GRASSMANN[pointlist,.95,5]

Making a small random perturbation of the boundary forces the algorithm to use the slit routine. This avoids the rounding error catastrophe but makes for slow progress for the algorithm.

pert = 
        Table[{Random[Real, {0,  .001}], Random[Real, {0,  .001}]}, {j, 1, 
            1000}];

points = pointlist +  pert;

GRASSMANN[points,.95,5]

A  snowflake

\(To\ demonstrate\ the\ generality\ of\ the\ Riemann\ mapping\ theorem\ it\ may\ be\ helpful\ to\ have\ an\ example\ with\ a\ more\ complicated\ boundary . \ 
          The\ Von\ Koch\ snowflake\ is\ approximated\ here\ out\ to\ five\ iterations\ in\ the\ construction . \ GRASSMAN\ performs\ reasonably\ well, \ 
      but\ some\ playing\ around\ with\ the\ example\ shows\ that\ the\ function\ "slitq"\ is\ \((\ again\ )\)\ in\ need\ of\ refinement . \ 
          For\ a\ much\ more\ \ efficient\ attack\ on\ the\ Von\ Koch\ snowflake\ see\ Don\ Marshall' s\ program\ ZIPPER\ \ at\) 
http: // www . math . washington . 
                edu/~marshall/zipper . html.

Clear[list]

Here are three points on the unit circle in an equillateral triangle.

list[
        0] = {{1, 0}, {- .5, Sqrt[3]/2}, {- .5, -Sqrt[3]/2}, {1, 
          0}}

The following function inserts the next three points in the construction between any two consecutive points from the prvious step.

step[pt1_, pt2_] := 
      N[{ (2/3) pt1  +  (1/3) pt2, (pt1  +  pt2)/
              2  -  (Sqrt[3]/6) {pt1[[2]] - pt2[[2]], 
                pt2[[1]] - pt1[[1]]},  (1/3) pt1  +  (2/
                  3) pt2, pt2}]

Now the kth iteration is defined recursively. For GRASSMANN to run, no points can be repeated in the list.

list[k_] := list[k] = 
        Join[{list[k - 1][[1]]}, 
          Flatten[Table[
              step[list[k - 1][[j]], list[
                    k - 1][[j + 1]]], {j, 1, Length[list[k - 1]] - 1}], 
            1]]

GRASSMAN will also have a problem if the points are spaced too far appart, so we fill in the polygonal figure with the vertices produced above.

fill[p1_, p2_] := 
      Table[(1 - t) p1 +  t p2, {t,  .01, 1,  .01}]

filledlist[list_] := 
      Join[{list[[1]]}, 
        Flatten[Table[
            fill[list[[j]], list[[j + 1]]], {j, 1, 
              Length[list] - 1}], 1]]

pointlist = filledlist[list[4]];

ListPlot[pointlist, AspectRatio -> Automatic]

The run time is a bit longer than for the previous examples.

GRASSMANN[pointlist,  .95, 1]

Conclusion

Grassman's algorithm probably has practical applications as a precursor for other more sophisticated methods .  The version presented here is obviously not meant for anything more than classroom demonstrations but could be altered to work on doubly connected regions or to work as a precursor  to produce input for the methods of Theodorsen or Wegmann  .  See [5]. The latter project would require the inverse correspondence produced in each step and so would present the usual numerical difficulties .  Readers interested in writing a better version  must of course refer to  Grassmann's paper [3] .   A footnote to [3] mentions that Grassmann was killed in a climbing accident in the year of publication .  It is satisfying to note that there appears to be not a single error in his  terse but detailed description and one would imagine that his implementation had many  fewer bugs than mine.

References

[1] L. Ahlfors, Complex  Analysis. An introduction to the theory of analytic functions of one complex variable. Third Edition.          
      McGraw-Hill Book Co., New York, 1978.

[2] Bishop, E., Foundations of constructive analysis. McGraw-Hill Book Co., New York-Toronto, Ont.-London, 1967

[3]  E.Grassmann, Numerical Experiments with a method of successive approximation for conformal mapping.
       Z. Angew .Math. Phys., 30, no.6, 1979, 873-884.

[4] P.Henrici, A general theory of osculation algorithms for conformal mapping. J. Lin. Alg. Appl. 52/53, 1983, 361-382

[5] P. Henrici, Applied and Computational Complex Analysis. Vol.3. John Wiley and Sons, Inc., New York, 1986.

[6] P. Koebe, Über eiene neue Methode der konformen Abbildung und Uniformisierung.,  Nachr. Kgl. Ges. Wiss. Göttingen, Math.-
      Phys. Kl. 1912, 844-848

[7 ] P. Koebe, Abhandlungen zur Theorie der konformen Abbildung,  J. Reine Angew. Math., 145, 1915, 177-223.

[8]  W. Osgood,  On the existence of Green's function for the most general simply connected plane region., Trans. Amer. Math.   
       Soc., 1 , 1900, 310-314 ,    http://www.jstor.org/view/00029947/di962766/96p0016u/

[9]  A. Ostrowski, Mathematische Miszellen XV: Zur konformen Abbildung einfach zusammenhängender Gebiete.
       Jahresber. Dt. Math.-Ver., 29, 1929, 168-182.

[10] R. Remmert, From Riemann surfaces to complex spaces. 
        Matériaux pour l'histoire des mathematiques au (XX^e) siécle, Nice, 1996, 203-241
        Sémin. Congr., 3, Soc. Math. France, Paris, 1998

[11] B. Riemann, Grundlagen für eine allgemeine Theorie der Funktionen einer veräderlichen complexen Grösse.
        Inauguraldissertation, Göttingen, 1851

[12]  P. Ullrich, The Riemann mapping problem. Studies in the history of Modern Mathematics, II. Rend. Circ. Mat. Palermo (2) 
        Suppl. No. 44, 1996, 9-42.

[13]  J.L.Walsh, History of the Riemann Mapping Theorem, Amer.Math.Monthly, 80,  no.3, 1973,  270-276
         http://www.jstor.org/view/00029890/di991580/99p0083d/0