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

Let the edges of an aperiodic Penrose tiling made from darts and kites be wires with finite resistances. By selecting outer and inner vertices as the current source and drains, calculate and visualize the potential and current distribution through the edge network. Use identical and random resistance values along the edges.

Off[General::"spell"]; Off[General::"spell1"];

(* the inflation rule *)
kiteDartRule = {p_, kd:("kite" | "dart")} :>
Module[{𝒻, a, b, c, d, e, f, g, h, i}, {a, b, c, d} = p;
       𝒻[a_, b_] = (a + GoldenRatio b)/(1 + GoldenRatio);
       (* construct new dart and kites *)
       Sequence @@ If[kd === "kite",
          e = 𝒻[a, c]; f = 𝒻[b, a]; g = 𝒻[d, a];
          h = a + c - d; i = a + c - b;
          {{{b, c, e, f}, "kite"}, {{d, g, e, c}, "kite"},
           {{f, e, a, h}, "dart"}, {{g, i, a, e}, "dart"}},
          (* last tile was a dart *) e = 𝒻[c, b]; f = 𝒻[c, d];
          g = b + c - a; h = d + c - a;
          {{{c, f, a, e}, "kite"}, {{e, a, b, g}, "dart"},
           {{f, h, d, a}, "dart"}}]];

(* start kite *)
pTiling[0] = {{{{0, 0}, GoldenRatio {Cos[36°], Sin[36°]},
               GoldenRatio {1, 0}, GoldenRatio {Cos[36°], -Sin[36°]}},
             "kite"}} // N[#, 30]&;
             
(* recursively apply kiteDartRule *)
pTiling[i_] := pTiling[i] = #[[1, 2]]& /@ 
 Split[Sort[{N[#[[1]]], #}& /@ (pTiling[i - 1] /. kiteDartRule)], 
            (* identical tiles? *) #1[[1]] === #2[[1]]&]

The functions edges and vertices give the edges and vertices after the oth iteration.

edges[o_] := edges[o] = Union[Sort /@ 
     Flatten[Partition[Append[#, First[#]], 2, 1]& /@ 
                         N[First /@ pTiling[o]], 1]];

vertices[o_] := vertices[o] = Union[Level[edges[o], {-2}]];

For indexing the voltages and currents, we prefer to use integers instead of approximate coordinates. So, we define the two function vertexNumber and vertexCoordinates that convert the explicit coordinates to integer numbers, back and forth.

makeVertexDefinitions[o_] :=
Do[(* cache coordinates and number of the vertices *)
   vertexNumber[o][vertices[o][[k]]] = k;
   vertexCoordinates[o][k] = vertices[o][[k]],
   {k, Length[vertices[o]]}]

allNumberedEdges[o_] := allNumberedEdges[o] =
           Map[vertexNumber[o], edges[o], {-2}];

Here are the counts of the vertices and edges for the first eight iterations. The eight iteration leads to more than 20000 equations for the current and voltages (edges and vertices).

Table[makeVertexDefinitions[o];
      {o, Length[edges[o]], Length[vertices[o]]}, {o, 8}] // 
 TableForm[#, TableHeadings -> {None, 
  StyleForm[#, FontWeight -> "Bold"]& /@ {"", "edges", "vertices"}}]&

Now we start with the first actual constellation. We use the 8th iteration and place the current source and drains at the left lower and right upper vertices.

order = 8;
{zeroVVertices[order], oneVVertices[order]} = List /@
vertexNumber[order] /@ {#[[1]], #[[-1]]}&[Sort[vertices[order]]]

Here is a sketch of this configuration.

Show[Graphics[{{Thickness[0.001], (* all edges *)
               {GrayLevel[1/2], Line[#]}& /@ edges[order]},
               (* source and drain vertices *)
               {PointSize[0.04], Hue[0], Point /@ 
                vertexCoordinates[order] /@ zeroVVertices[order]},
               {PointSize[0.04], Hue[0.76], Point /@ 
                vertexCoordinates[order] /@ oneVVertices[order]}}],
                  AspectRatio -> Automatic]

To write down the current conservation law at all vertices (the source and drain vertices are held at a fixed potential), we construct list of all edges that meet at a given vertex.

edgesAtVertices[o_] := edgesAtVertices[o] =
{#[[1, 1]], Last /@ #}& /@
Split[Sort[Flatten[{{#[[1]], #}, {#[[2]], #}}& /@ 
           allNumberedEdges[o], 1]], #1[[1]] === #2[[1]]&];

Now, we write all equations governing the current flow through the edges of the tiling. Along all edges, we have i[a, b] = (u[b] - u[a])/R[a, b] where i[a, b] is the current from vertex a to vertex b. And with all edges oriented, we can write down the current conservation equations for all nonsource and nondrain vertices. For a given network specified by the edges and the connectivity, the function makeUIEquations returns the equations for the voltage values at the vertices and for the currents through the edges. Without loss of generality (because we assume a linear relation between current and voltage), we set the source potential to zero and the drain potential to one.

makeUIEquations[numberedEdges_, edgesAtVertices_,
                {zeroVVertices_, oneVVertices_}, R_, {u_, i_}] :=
Flatten @ {
(* voltage drops across a resistor *)
i[#1, #2] == (u[#2] - u[#1])/R[#1, #2]& @@@ 
          numberedEdges,
(* current conservation at a vertex *)
Function[{vertex, edgeList}, 0 == (* all directed currents *)
 Plus @@ ((If[# === vertex, 1, -1] i[##])& @@@ edgeList)] @@@ 
    DeleteCases[edgesAtVertices, 
                (* no current conservation for 
                   source and drain vertices *)
                Alternatives @@ ({#, _}& /@ 
                Flatten[{zeroVVertices, oneVVertices}])]} /.
    (* source and drain voltage values *)
    Dispatch[Join[(u[#] -> 0)& /@ zeroVVertices, 
                  (u[#] -> 1)& /@ oneVVertices]]

We start by using identical resistance (value 1) for all edges. Here are the 20262 equations in 20262 variables.

eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                      {zeroVVertices[order], oneVVertices[order]}, 
                      1&, {u, i}];

iuVariables[eqs_] := Cases[eqs, _u | _i, Infinity] // Union

vars = iuVariables[eqs];
Length /@ {eqs, vars}

We solve the equations with a direct call to Solve. For an efficient solution, we numericalize the equations.

(* avoid numericalizations of the variable numbering *)
SetAttributes[{u, i}, NHoldAll];
sol1 = Solve[N @ eqs, vars];

To visualize the result sol1, we implement the two functions voltagePlot and currentPlot. The first plot draws the vertices colored according to the voltage. Blue is low and red is high. The second is a plot of the magnitudes of the current along the edges. Again blue is lowest and red is highest.

voltagePlot[o_, sol_, u_] :=
Module[{λ = Cases[sol, _u, Infinity] // Union // Length, us, cs},
(* different voltage values *)
us = Union[Rationalize[SetPrecision[#, 6], 0]& @
     Sort[Last /@ Cases[sol, HoldPattern[_u -> _], Infinity]]];
(* make groups of voltage values and color function *)
cs = Round[Length[us]/10]; ipo = Interpolation[
MapIndexed[{#1, #2[[1]] - 1}&, Flatten[{0, 
    #[[Round[Length[#]/2]]]& /@ (* homogeneous coloring *)
                Partition[us, Round[Length[us]/cs]], 1}]]];
(* make graphics *)
Show[Graphics[{(* scale points *) PointSize[1/Sqrt[λ]/2], 
({Hue[0.78 (1 -  Abs[ipo[u[#]]/cs])], 
Point[vertexCoordinates[o][#]]}& /@  
      Range[Length[vertices[o]]]) /. Dispatch[sol[[1]]] /.
    (* the drain and source vertices *)
    Dispatch[Join[(u[#] -> 0)& /@ zeroVVertices[o], 
                  (u[#] -> 1)& /@ oneVVertices[o]]]}],
      AspectRatio -> Automatic]];

currentPlot[o_, sol_, i_] :=
Module[{λ = Length[Union[Cases[sol, _i, Infinity]]], is, cs},
(* different current values *)
is = Union[Rationalize[SetPrecision[#, 6], 0]& @
     (is0 = Sort[Abs[Last /@ Cases[sol, HoldPattern[_i -> _],
                 Infinity]]])];
(* make groups of voltage values and color function *)
cs = Round[Length[is]/20]; ipo = Interpolation[
MapIndexed[{#1, #2[[1]] - 1}&, Flatten[{is0[[1]], 
#[[Round[Length[#]/2]]]& /@ 
 Partition[is, Round[Length[is]/cs]], is0[[-1]]}]]];
(* make graphics *)
Show[Graphics[{(* scale line thickness *) Thickness[1/λ], 
({Hue[0.78 (1 - Abs[ipo[Abs[i @@ #]]]/cs)], 
  Line[vertexCoordinates[o] /@ #]}& /@ allNumberedEdges[o]) /. 
                                  Dispatch[sol[[1]]]}],
      AspectRatio -> Automatic]]

We unite the two last types of plots in the function voltageCurrentPlot.

(* show voltages and currents together *)
voltageCurrentPlot[o_, sol_] := 
Show[GraphicsArray[
Block[{$DisplayFunction = Identity},
      {voltagePlot[o, sol, u], currentPlot[o, sol, i]}]]];

Here are the resulting voltage and current distribution for our first example. Coarse grained, the voltage is an increasing function from the left lower to the right upper corner. (The voltages can be considered as a discrete approximation to the Laplace equation of the electric potential. And, as a harmonic function it cannot have a local minimum or maximum.) The currents flow mainly along the diagonal. The occasional blue edges along the diagonal area are due to nearly or exactly identical voltages at their endpoints.

voltageCurrentPlot[order, sol1];

We repeat the voltage and current calculation for the last network with random resistances. We start with a uniform distribution of the resistance values. Now, the current distribution looks less regular.

SeedRandom[1];
(* form equations *)
eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                     {zeroVVertices[order], oneVVertices[order]}, 
                     Random[]&, {u, i}];
(* visualize solutions *)                     
voltageCurrentPlot[order, sol2 = Solve[N @ eqs, vars]]

Next, we use a more stretched out distribution with resistance values spanning eight orders of magnitude. The current from the left lower vertex to the right upper vertex occurs now along more pronounced paths instead in a more uniform manner along the diagonal. The voltage distribution looks nearly identical to the last one.

SeedRandom[2];
(* form equations *)
eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                     {zeroVVertices[order], oneVVertices[order]}, 
                     10^Random[Real, {-4, 4}]&, {u, i}];
(* visualize solutions *)                      
voltageCurrentPlot[order, sol3 = Solve[N @ eqs, vars]]

At this point, let us check the quality of the solutions. The maximal residual is about 2×(10^(-11)).

Max[Abs[(Subtract @@@ eqs) /. Dispatch[sol3[[1]]]]]

We could now continue by making the resistance of some edges infinite by using If[Random[] < p, 1, Infinity], effectively reducing the number of potential current paths from the source to the drain, but leave this to the reader.

We continue with a configuration where we make all outer vertices are current drains and the central vertex is the source. The functions outerVertices and centralVertex calculate these vertices for a given iteration order o.

outerVertices[o_] := outerVertices[o] = 
vertexNumber[order] /@ Union[Flatten[
Select[Split[Sort[Sort /@ Flatten[
              Partition[Append[#, First[#]], 2, 1]& /@ 
                       N[First /@ pTiling[o]], 1]]],
       Length[#] == 1&], 2]];

centralVertex[o_] := centralVertex[o] =
Module[{mp = Plus @@ vertices[o]/Length[vertices[o]], distances},
       distances = Sqrt[#.#]&[mp - #]& /@ vertices[o];
       vertexNumber[o] @ vertices[o][[
              Position[distances, Min[distances]][[1, 1]]]]]

(* the outer and inner vertices for the current iteration order *)
{zeroVVertices[order], oneVVertices[order]} = 
{outerVertices[order], {centralVertex[order]}};

Here are the resulting voltage and current distributions for all edges having identical resistance values.

(* form equations *)
eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                     {outerVertices[order], {centralVertex[order]}}, 
                     1&, {u, i}];

vars = iuVariables[eqs];
(* visualize solutions *) 
voltageCurrentPlot[order, sol4 = Solve[N @ eqs, vars]]

Using again random resistance values along the edges shows a network of current paths. The voltage distribution looks similar but has some different details.

(* use two different resistance values *)
SeedRandom[3];
(* form equations *)
eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                     {outerVertices[order], {centralVertex[order]}}, 
                      If[Random[] < 0.1, 1, 10^6]&, {u, i}];
(* visualize solutions *) 
voltageCurrentPlot[order, sol5 = Solve[N @ eqs, vars]]

(* use continuous distribution of resistance values *)
SeedRandom[4];
(* form equations *)
eqs = makeUIEquations[allNumberedEdges[order], edgesAtVertices[order],
                     {outerVertices[order], {centralVertex[order]}}, 
                      10^Random[Real, {-4, 4}]&, {u, i}];
(* visualize solutions *) 
voltageCurrentPlot[order, sol6 = Solve[N @ eqs, vars]]

We end with a plot of the currents on a logarithmic scale.

Show[
Graphics[{Hue[#1], PointSize[0.002],
  MapIndexed[If[#1 == 0., {}, Point[{#2[[1]], Log[10, #1]}]]&,
      Sort[Abs[Last /@ Cases[#2, HoldPattern[_i -> _], 
     Infinity]]]]}]& @@@
   {{0, sol1}, {0.1, sol2}, {0.2, sol3}, {0.33, sol4}, 
    {0.44, sol5}, {0.78, sol6}},
Frame -> True, PlotRange -> {-15, 0}]

For some theoretical results about the current distribution of random resistor networks, see [✶], [✶], [✶], [✶], [✶]. For the resistance of a Sierpinski carpet, see [✶].