This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

The Schwarz–Christoffel map from the unit disk to a nonintersecting simple polygon 𝒫 with vertices (Subscript[w, k]) is given by [✶], [✶]

f (z)=α+β(∫_0^z)(∏_(k=1)^n)(((1-(ζ/(Subscript[z, k]))))^((Subscript[α, k])-1))dζ

where f(z)∈𝒫 for |z|≤1 and f((Subscript[z, k]))=(Subscript[w, k]) and (Subscript[α, k]) π are the counterclockwise-oriented interior angles of the polygon. Calculate a numerical approximation to a Schwarz–Christoffel map for an L-shaped polygonal domain with vertices {0,2,2+i,1+i,1+2i,2i}. Visualize the resulting conformal map.

wList = {0, 2, 2 + I, 1 + I, 1 + 2 I, 2 I} - (1 + I)/2;

αList = {Pi/2, Pi/2, Pi/2, 3Pi/2, Pi/2, Pi/2}/Pi;

Show[Graphics[
  {Thickness[0.02], Line[{Re[#], Im[#]}& /@ Append[wList, First[wList]]],
   {Hue[0], PointSize[0.04], Point[{0, 0}]}},
    PlotRange -> All, AspectRatio -> Automatic, Frame -> True]]

Calculating the Schwarz–Christoffel map boils down to determining the preimages (Subscript[z, k]) of the vertices and the two constants α and β. The last two constants can be determined by prescribing the position of two vertices. We will use the two vertices on the symmetry axis Re(z)=Im(z) (the first and the fourth vertex). And taking the symmetry of the L-shaped domain into account, the only two unknowns are the points (Subscript[z, 2]) and (Subscript[z, 3]). Because the boundary of the unit disk gets mapped to the boundary of the polygon, we write (Subscript[z, k])=exp(i (Subscript[φ, k])) and we are left with determining the two values (Subscript[φ, 2]) and (Subscript[φ, 3]). For given values of  (Subscript[φ, 2]) and (Subscript[φ, 3]), the function makeSymmetriczList gives the list of all six vertices of the L-shaped polygon.

makeSymmetriczList[{φ2_, φ3_}] = 
{Exp[I 5Pi/4], Exp[I φ2], Exp[I φ3], Exp[I Pi/4],
 Exp[I Pi/4] Conjugate[Exp[-I Pi/4] Exp[I φ3]],
 Exp[I Pi/4] Conjugate[Exp[-I Pi/4] Exp[I φ2]]};

Using numerical integration, the function wValues calculates the positions of the six vertices for given values (Subscript[φ, 2]) and (Subscript[φ, 3]).

wValues[{φ2_, φ3_}] :=
Module[{zs, ws, αβSol},
 zs = makeSymmetriczList[{φ2, φ3}];
 ws = Table[NIntegrate[Evaluate @ (Times @@ ((1 - ζ/zs)^(αList - 1))), 
                       {ζ, 0, zs[[j]]}], {j, 6}];
 (* position first and fourth vertex through α and β *)
 αβSol = Solve[{wList[[1]] == α + β ws[[1]], 
                wList[[4]] == α + β ws[[4]]}, {α, β}][[1]];
 (* return linear map and mapped points *)
 {αβSol, α + β ws /. αβSol}]

Minimizing now the sum of squares of the distances of the mapped preimages of the vertices should yield a minimum with value 0 and the values of (Subscript[φ, 2]) and (Subscript[φ, 3]) for the Schwarz–Christoffel map.

g[{φ2_?NumberQ, φ3_?NumberQ}] := Plus @@ (Abs[wList - wValues[{φ2, φ3}][[2]]]^2)

fm = FindMinimum[g[{φ2, φ3}], {φ2, -1, -1.1},  {φ3, 0.5, 0.4}]

So, we obtain the following parameters and integrand for the map.

(* constants α and β for the final map *)
{αF, βF} = {α, β} /. wValues[{φ2, φ3} /. fm[[2]]][[1]];

(* the final z_k values *)
finalZs = makeSymmetriczList[{φ2, φ3} /. fm[[2]]];

(* the final integrand *)
scM[ζ_] = N[(Times @@ ((1 - ζ/finalZs)^(αList - 1)))]

To visualize the map we will not use NIntegrate, but rather NDSolve. For a given direction 𝒹, we solve the differential equation corresponding to the integral in f(z) radially outward.

solveODERadial[𝒹_] := NDSolve[{w'[t] == 𝒹 scM[𝒹 t], w[0] == 0}, w, 
                              {t, 0, 1 - 10^-10}]

(* radial distances for equi-lines *)
pp = 24;
tPoints = Join[Table[t, {t, 0, 1 - 1/pp, 1/pp}], {0.986, 0.9999}];

data = Table[ndsol = solveODERadial[Exp[I φ]];
             ((αF + βF w[#])& /@ tPoints) /. ndsol[[1]], 
             {φ, 0, 2Pi, 2Pi/360}];

(* the mapped points *)
ε = 10^-6;
data = Table[ndsol = solveODERadial[Exp[I φ]];
             ((αF + βF w[#])& /@ tPoints) /. ndsol[[1]], 
             {φ, Pi/4 + ε, 9/4 Pi - ε, (2Pi - 2 ε)/360}];

Show[Graphics[
{(* the L-shaped dpolygon *)
 {Hue[0], Thickness[0.006], Line[{Re[#], Im[#]}& /@ 
  Append[wList, First[wList]]]},
 (* equi-lines and perpendicular curves *)
 {GrayLevel[0.5], MapIndexed[If[IntegerQ[#2[[1]]], 
                                Line[{Re[#], Im[#]}& /@ #], {}]&, data]},
 {GrayLevel[0], Line[{Re[#], Im[#]}& /@ #]& /@ Transpose[data]}}], 
      PlotRange -> All, Frame -> True, AspectRatio -> Automatic]

The last graphic shows the effect of crowding: near the left lower corner, the curves perpendicular to the equi-lines cluster and in the arms the become relatively widely spaced. To get a more homogeneously spaced equi-line graphic, we will now solve the differential equation azimuthally and use ParametricPlot to get adaptively refined smooth equi-curves.

solveODEAzimuthal[ρ_, {φ1_, φ2_}] := 
With[{w0 = αF + βF NIntegrate[scM[ζ], {ζ, 0, ρ Exp[I φ1]}]}, 
     NDSolve[{w'[φ] == βF ρ I Exp[I φ] scM[ρ Exp[I φ]], w[φ1] == w0}, w, 
             {φ, φ1, φ2}]]

contract[{a_, b_}] := With[{mp = (a + b)/2, lh = (b - a)/2, δ = 1 - 10^-10},
                           {mp - δ lh, mp + δ lh}]

φRanges = contract /@ 
           Partition[Append[#, 2Pi + First[#]]&[Arg[finalZs // N]], 2, 1];

The function SmoothPlotAzimuthal generates the plots from the list of  InterpolatingFunction returned by NDSolve.

makeSmoothPlotAzimuthal[ρ_, ipo_InterpolatingFunction] := 
ParametricPlot[{Re[ipo[φ]], Im[ipo[φ]]}, Evaluate[Flatten[{φ, ipo[[1, 1]]}]],
               PlotPoints -> 120, DisplayFunction -> Identity]

makeSmoothPlotAzimuthal[ρ_, l_List] := 
    Show[makeSmoothPlotAzimuthal[ρ, #]& /@ 
                   Cases[l, _InterpolatingFunction, Infinity]]

We also display some of the perpendicular curves, namely the ones ending at the vertices of the L-shaped polygon and at the midpoints of its edges.

makeSmoothPlotRadial[𝒹_, ipo_InterpolatingFunction] := 
ParametricPlot[{Re[αF + βF ipo[t]], Im[αF + βF ipo[t]]}, 
               Evaluate[Flatten[{t, ipo[[1, 1]]}]],
               PlotPoints -> 120, DisplayFunction -> Identity,
               PlotStyle -> {{GrayLevel[0.5], Thickness[0.003]}}]

makeSmoothPlotRadial[𝒹_, l_List] := Show[makeSmoothPlotRadial[𝒹, #]& /@ 
                             Cases[l, _InterpolatingFunction, Infinity]]

(* calculate the mapped point position *)
of[φ_?NumberQ] := Module[{res},
                         Off[NIntegrate::ncvb];
                         res = αF + βF NIntegrate[scM[ζ], {ζ, 0, Exp[I φ]}];
                         On[NIntegrate::ncvb]; res];

(* calculate z-values corresponding to the midpoints of the L edges *)
edgeMidpointZs = 
(Exp[I Re[φ]] /. FindRoot[of[φ] - #1, {φ, #2[[1]], #2[[2]]},
                          PrecisionGoal -> 8])& @@@ 
Transpose[{(Plus @@@ Partition[Append[wList, First[wList]], 2, 1]/2.),
           Partition[Append[#, 2Pi + First[#]]&[Arg[finalZs // N]], 2, 1]}]

Along the curve to the midpoint of the second polygon edge, we calculate equi-line radii that are approximately homogeneously spaced.

𝒻[t_] = αF + βF (w /. solveODERadial[edgeMidpointZs[[2]]][[1]])[t];

i𝒻 = Interpolation[Table[{Re[𝒻[t]], t}, {t, 0, 1, 1/12}]];

ρs = Table[i𝒻[t], {t, Re[𝒻[0]], Re[𝒻[1]], (Re[𝒻[1]] - Re[𝒻[0]])/15}]

Now we obtain smooth, approximately homogeneously spaced equicurves up to near the polygon boundary.

Show[{Graphics[{Hue[0], Thickness[0.01], Line[{Re[#], Im[#]}& /@ 
                Append[wList, First[wList]]]}], 
      (* radial orthogonal curves *)
      Function[𝒹, makeSmoothPlotRadial[𝒹, solveODERadial[𝒹]]] /@ 
                                   Join[finalZs, edgeMidpointZs],
      (* smooth equi-cuves *)
      Function[ρ, makeSmoothPlotAzimuthal[ρ, 
                  solveODEAzimuthal[ρ, #]& /@ φRanges]] /@ 
                  (* radii for homogeneously spaced equi-lines *)           
                                         Append[Most[ρs], 0.998]},
          PlotRange -> All, AspectRatio -> Automatic,
          DisplayFunction -> $DisplayFunction];

We end with a visualization of the Schwarz–Christoffel map of a more complicated polygon.

(* vertices of a polygon in the complex plane *)
wList = With[{s = Sqrt[3]}, {1 + I s, 2 + 2 I s, 1 + 3 I s, -3 + 3 I s, 
               -2 + 2 I s, 2 I s, -2, -4, -5 - I s, -3 - 3 I s, -2 - 2 I s, 
               -3 - I s, 1 - I s, 2 - 2 I s, 4 - 2 I s, 6, 4, 3 - I s}/2];

equiAndFieldCurves[equiValues_List] := 
Module[{αList, φs, φEs, φL, β, zs, w0, eqs, ρ, ϕ1, ϕ2, ndsol, pp, 
        curve, ℛ, ω, opts, R = 1 - 10^-12},
 (* angles in units of π of the polygon *)
 αList =  1/3 {4, 2, 2, 1, 2, 5, 4, 2, 2, 1, 2, 5, 4, 2, 2, 1, 2, 5};       
 (* parameters for the Schwarz-Christoffel map *)
 φs = {1.880275011280091, 2.074403517501587, 2.084559146579438, 
       2.087831881871811, 2.087952818861801, 2.094395102393195};
 (* parameters for edge midpoints *)
 φEs = {2.056709598824023, 2.0812111099343658, 2.0875447164599064,
        2.087856604171078, 2.0882552129179253, 2.5163592862103163};
 β = 0.813037054574698 - 0.56337234858643208 I;
 (* rotation matrices *)
 Do[ℛ[k] = N[{{Cos[k 2Pi/3], Sin[k 2Pi/3]}, {-Sin[k 2Pi/3], Cos[k 2Pi/3]}}], 
    {k, 0, 2}];
 (* graphics options; allow unlimited subdivision *)
 optsG = Sequence[PlotPoints -> 100, PlotDivision -> Infinity, 
                     MaxBend -> 3, DisplayFunction -> Identity];
 (* NDSolve options *)
 optsD = Sequence[PrecisionGoal -> 12, AccuracyGoal -> 12, 
                  WorkingPrecision -> 25];                    
 (* all preimages of the 18 vertices *)
 zs = Join[Exp[I φs], Exp[2Pi I/3] Exp[I φs], Exp[4Pi I/3] Exp[I φs]];
 (* equi- and filed lines *)
 {(* equilines *)
  Table[ρ = equiValues[[j]]; ϕ1 = φs[[1]]; ϕ2 = ϕ1 + 2/3 Pi;
  (* starting value for differential equations *)
  w0 = β NIntegrate[Evaluate @ SetPrecision[(Times @@ ((1 - ζ/zs)^(αList - 1))), 30], 
                    Evaluate @ {ζ, 0, SetPrecision[ρ zs[[1]], 30]},
                    MaxRecursion -> 20, PrecisionGoal -> 14, WorkingPrecision -> 30];
  (* solve differential equation azimuthally *)
  eqs = {1/(I ρ Exp[I φ]) Y'[φ] == β (Times @@ ((1 - ρ Exp[I φ]/zs)^(αList - 1))),
         Y[ϕ1] == w0};
  ndsol = NDSolve[SetPrecision[eqs, 30], Y, {φ, ϕ1, ϕ2}, optsD];
  (* make plot of curve *)
  pp = ParametricPlot[Evaluate[{Re[Y[φ]], Im[Y[φ]]} /. ndsol[[1]]], 
                     {φ, ϕ1, ϕ2}, Evaluate[optsG]];
  (* rotate curve in other positions *)
  curve = Cases[pp, _Line, Infinity];
  Table[Map[ℛ[k].#&, curve, {-2}], {k, 0, 2}], {j, Length[equiValues]}],
 (* field lines from the origin to the 18 polygon vertices *) 
 φL = Join[φs, φEs];
 Table[ω = φL[[j]]; Clear[ρ];
 eqs = {1/(Exp[I ω]) Y'[ρ] == β (Times @@ ((1 - ρ Exp[I ω]/zs)^(αList - 1))),
        Y[0] == 0};
  (* solve differential equation radially *)
  ndsol = NDSolve[SetPrecision[eqs, 30], Y, {ρ, 0, R}, optsD];
  (* make plot of curve *)
  pp = ParametricPlot[Evaluate[{Re[Y[ρ]], Im[Y[ρ]]} /. ndsol[[1]]], 
                     {ρ, 0, R}, Evaluate[optsG]];
  curve = Cases[pp, _Line, Infinity];
  (* rotate curve in other positions *)
  Table[Map[ℛ[k].#&, curve, {-2}], {k, 0, 2}], {j, 12}]}]

Module[{(* equi-values to be used *)
        equiList = {0.1, 0.3, 0.5, 0.7, 0.88, 0.96, 0.99, 0.997, 
                    0.999, 0.9996, 0.99992, 0.99997, 0.9999988}, λ, 
                    equiCurves, fieldCurves},
       λ = Length[equiList];
       (* equi- and field curves *)
       {equiCurves, fieldCurves} = equiAndFieldCurves[equiList];
       (* make graphics *)
       Show[Graphics[
        {(* polygon in thick gray *)
         {Thickness[0.008], GrayLevel[1/2], 
          Line[{Re[#], Im[#]}& /@ Append[wList, First[wList]]]}, 
         (* field curves in black *)
         {Thickness[0.002], GrayLevel[0], fieldCurves},
         (* colored equi-curves *)
         {Thickness[0.002], MapIndexed[{Hue[0.78 (#2[[1]] - 1)/(λ - 1)], #}&, 
                                        equiCurves]}}], 
             AspectRatio -> Automatic, PlotRange -> All]];