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

Take the 2D Hilbert curve from Exercise 20 of Chapter 2 of the Graphics volume and embed it in a square grid (so that the edges of the Hilbert curve are along neighboring lattice points). Assuming springs of stiffness 1 between all neighboring grid points and springs of stiffness 𝓀 along the segments of the Hilbert curve, allow the point mass lattice points to oscillate in vertical direction (outermost lattice points are kept fixed). Calculate and visualize some of the resulting eigenoscillations for various stiffnesses 𝓀.

LSystemWithFAndLAndR[axiom_, rules_, iterations_, δ_] :=
Module[{minus, plus, ℛ, last = {0, 0}, d = {1, 0}},
 minus = {{ Cos[ δ], Sin[ δ]}, {-Sin[ δ], Cos[ δ]}};
 plus  = {{ Cos[-δ], Sin[-δ]}, {-Sin[-δ], Cos[-δ]}};
 ℛ = rules /. {(a_ -> b_) -> (a :> Sequence @@ b)};
Select[Prepend[
 (Which[# == "F", last = last + d, # == "+", d = plus.d; Null,
        # == "-", d = minus.d; Null, # == "R", Null,   
        # == "L", Null]& /@ Nest[(# /. ℛ)&, axiom, iterations]),
      {0, 0}], # =!= Null&]]

order = 4; {L =  3^(order) + 1, dim = 3^(2 order)}

Show[Graphics[
 MapIndexed[{Hue[(#2[[1]] - 1)/3^(2 order)],
  Thickness[0.001], Line[#]}&, Partition[line = 
  {1, 1} + {1, -1}#& /@ LSystemWithFAndLAndR[{"-", "L"},
   {"L" -> Characters["LF+RFR+FL-F-LFLFL-FRFR+"],
    "R" -> Characters["-LFLF+RFRFR+F+RF-LFL-FR"]},
    order, 2Pi/4], 2, 1]]], AspectRatio -> Automatic, 
     Frame -> True, PlotRange -> L {{0, 1}, {0, 1}}]

To calculate the oscillations of the vertical spring network, we need the connectivity of the grid. The equation of motion for the lattice points (Subscript[u, i,j])(t) is (Subscript[(Overscript[u, ..]), i,j])(t)=∑(Subscript[k, i,j]^(i',j'))((Subscript[u, i',j'])(t)-(Subscript[u, i,j])(t)) where (Subscript[k, i,j]^(i',j')) is the stiffness of the spring that connects the lattice points {i,j} and {i',j'} and the sum extends over the four neighboring lattice points of {i,j} . Assuming (Subscript[(Overscript[u, ..]), i,j])(t)∼(e^(i ω t))(Subscript[(OverHat[u]), i,j]), the calculation of the eigenoscillations reduces to the calculation of the eigenvalues and eigenvectors of the coefficient matrix of ∑(Subscript[k, i,j]^(i',j'))((Subscript[(OverHat[u]), i',j'])-(Subscript[(OverHat[u]), i,j])).

To obtain a linear numbering of the lattice points, we number them consecutively, first in the x-, then in the y-direction.

Do[vertexNumber[{i, j}] = (j - 1) (L - 1) + i, {i, L -1}, {j, L - 1}]

Next, we find all lattice points that are connected by the Hilbert curve.

hData = Function[x, {x[[1, 1, 1]], {#1[[2]], #2}& @@@ x}] /@
Split[Sort[Flatten[{{#1, #2}, {Reverse[#1], #2}}&  @@@
 MapIndexed[{#1, 𝓀[#2[[1]]]}&, Partition[line, 2, 1]], 1]],
      (#[[1, 1]] === #2[[1, 1]])&];

Now, we assemble the coefficient matrix. It is a sparse 6561×6561 matrix with 32481 elements.

(* make rules for coefficients for neighbors *)
makeRules[{p_, ps_List}] :=
Block[{qs, pn = vertexNumber[p], ℓ = Length[ps]}, 
 qs = DeleteCases[
      Complement[(p + #)& /@ {{-1, 0}, {1, 0}, {0, -1}, {0, 1}},
                 ps], {0, _} | {L, _} | {_, 0} | {_, L}];  
 (* the nonvanishing matrix elements *)
 Flatten[{{pn, pn} -> (Plus @@ (Last /@ ps)) + (4 - ℓ) 1,
          ({pn, vertexNumber[#1]} -> -#2)& @@@ ps,
          ({pn, vertexNumber[#1]} -> -1)& /@ qs}]]

(* all coefficients *)
allRules[𝓀_] = makeRules /@ hData;
(* create sparse arrays *)
M[𝓀_] := SparseArray[Flatten[allRules[𝓀]]]

We use three different values of 𝓀: a small one, a large one, and 𝓀=1 (meaning a uniform spring network). Here are the corresponding values of extreme eigenvalues.

(* three stiffness values *)
𝓀List = {1/100, 1, 100};
Evaluate[ℳ /@ 𝓀List] = N[M /@ (Function /@ 𝓀List)];

(εRange[#] = 
Reverse[Eigenvalues[ℳ[#], 2, 
  Method -> {"Arnoldi", Criteria -> "BothEnds",  BasisSize -> 20,
             MaxIterations -> 10^4}]])& /@ 𝓀List

Given a list of eigenvectors, the function makeGraphicsArrays creates graphics arrays with the eigenoscillations of the three different spring stiffnesses adjacent.

makeGraphicsArrays[evsList_] := Show[GraphicsArray[#]]& /@ 
Map[ListDensityPlot[Partition[#, L - 1], Mesh -> False, 
       PlotRange -> All, DisplayFunction -> Identity,
       FrameTicks -> None]&, Transpose[evsList], {-2}]

Finally, we set some option values for the Arnoldi iterations (a smaller than the default tolerance is enough for visualization purposes). We use the Shift suboption to obtain eigenfunctions to different elements of the spectrum.

𝒪[λ___] := Method -> {"Arnoldi", BasisSize -> 40, 
              MaxIterations -> 10^4, Tolerance -> 10^-5,
              If[{λ} =!= {}, Shift -> λ, Sequence @@ {}]}

Here are the lowest states for the three spring networks. On the background of the 𝓀=1 springs, we see modulations due to the different spring stiffnesses along the Hilbert curve.

lowestStates = Reverse[Eigenvectors[ℳ[#], -4, 𝒪[]]& /@ 𝓀List];
makeGraphicsArrays[lowestStates];

The highest states oscillate as frequently as possible.

highestStates = Eigenvectors[ℳ[#], 4, 𝒪[]]& /@ 𝓀List;
makeGraphicsArrays[highestStates];

The eigenfunctions corresponding to central eigenvalues show a complicated behavior. The imprint of the Hilbert curve springs is clearly visible.

Do[Print["α = ", α];
middleStates = Eigenvectors[ℳ[#], 3, 𝒪[#2]]& @@@ Transpose[
                  {𝓀List, (#1 + α (#2 - #1))& @@@ (εRange /@ 𝓀List)}];
makeGraphicsArrays[middleStates], {α, 1/4, 3/4, 1/4}]

We end with two random realizations where each individual spring along the Hilbert curve has a random stiffness.

SeedRandom[1];
{ℳR[1], ℳR[2]} = (* realize random spring constant *)
Table[Block[{𝓀},  𝓀[k_] := 𝓀[k] = 10^Random[Real, {-1, 1}]; M[𝓀]], {2}];

(* lowest and highest eigenvalues *)
(εRange[#] = 
Reverse[Eigenvalues[ℳR[#], 2, 
  Method -> {"Arnoldi", Criteria -> "BothEnds",  BasisSize -> 20,
             MaxIterations -> 10^4}]])& /@ {1, 2}

The lowest states are not changed globally, but the higher states are become quickly localized due to the random spring stiffnesses.

(* lowest states *)
makeGraphicsArrays[Reverse @ Eigenvectors[ℳR[#], -2, 𝒪[]]& /@ {1, 2}]

(* intermediate states *)
makeGraphicsArrays[Reverse @ Eigenvectors[ℳR[#], +4, 𝒪[4]]& /@ {1, 2}]

(* highest states *)
makeGraphicsArrays[Reverse @ Eigenvectors[ℳR[#], +4, 𝒪[8]]& /@ {1, 2}]