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.

Symbolics Volume, Chapter 1

Consider a level three Gosper curve from Chapter 2 of the Graphics volume. Assuming harmonically bound point masses at the vertices (the two outermost points are kept fixed) and springs of a given stiffness along the straight line segments of the Gosper curve, the resulting object is a 2D vibrating system. For small elongations, calculate the resulting eigenoscillations for identical stiffnesses and random stiffnesses. Starting with one inner vertex initially elongated, visualize how the vibrations propagate through the spring chain in the two cases.

LSystemWithlAndr[axiom_, rules_, iterations_, δ_] :=
Module[{m, p, c = Cos[δ], s = Sin[δ], ℛ, ℓ = {0, 0}, 𝒹 = {1, 0}},
 m = {{ c, s}, {-s, c}}; p = {{ c, -s}, { s, c}};
 ℛ = rules /. {(a_ -> b_) -> (a :> Sequence @@ b)};
 Select[Prepend[
 (Which[# == "l", ℓ = ℓ + 𝒹,  # == "r", ℓ = ℓ + 𝒹,
        # == "+", 𝒹 = p.𝒹; Null, # == "-", 𝒹 = m.𝒹; ]& /@
       Nest[(# /. ℛ)&, axiom, iterations]), {0, 0}], # =!= Null&]]

Gosper[i_Integer] := Gosper[i] = LSystemWithlAndr[
{"l"}, {"l" -> Characters["l+r++r-l--ll-r+"],
        "r" -> Characters["-l+rr++r+l--l-r"]}, i, Pi/3]

Next, we join the adjacent parallel line segments into one.

(* are the three points p1, p2, p3 on one straight line? *)
oneOneLineQ[{p1_, p2_, p3_}] := oneOneLineQ[{p1, p2, p3}] =
With[{𝓃 = Sqrt[#.#]&}, 
     Expand[𝓃[p3 - p1] - (𝓃[p2 - p1] + 𝓃[p3 - p2])] === 0]

order = 3;
gosperLine = Gosper[order] //. {a___, p1_, p2_, p3_, b___} :>
              {a, p1, p3, b} /; oneOneLineQ[{p1, p2, p3}];

Here is the resulting Gosper curve with 293 inner vertices. We color each line segment differently. For later use, we mark one central point.

{𝒹 = Length[gosperLine], dim = 𝒹 - 2}

Position[#, Min[#]]&[
Function[mp, (mp - #).(mp - #)& /@ gosperLine][
         Plus @@ gosperLine/Length[gosperLine]]][[1, 1]]

(* position of one central vertex *)
m = Which[order == 1, 4, order == 2, 15, 
          order == 3, 113, order == 4, 791];
(* or use point nearest to center:
Position[#, Min[#]]&[
Function[mp, (mp - #).(mp - #)& /@ gosperLine][
         Plus @@ gosperLine/Length[gosperLine]]][[1, 1]] *)          
(* show straight line segments of the Gosper curve *)
Show[Graphics[{{Hue[Random[]], Thickness[0.002], Line[#]}& /@ 
                Partition[gosperLine, 2, 1],
               {PointSize[0.03], Point[gosperLine[[m]]]}}], 
     AspectRatio -> Automatic]

Giving the first vertex the index 0, allows to conveniently refer to the 293 movable vertices by their number vertexNumber.

Do[vertexNumber[gosperLine[[k]]] = k - 1;
   vertexCoordinates[k - 1] = gosperLine[[k]],
   {k, Length[gosperLine]}]

To calculate the eigenoscillations of the Gosper spring chain, we make the ansatz (Subscript[(Overscript[u, ..]), j])(t)∼(e^(i ω t))(Subscript[(OverHat[u]), j]) in the equations of motions (Subscript[(Overscript[u, ..]), j])(t)=(Subscript[k, j ])f((Subscript[u, j-1])(t),(Subscript[u, j])(t))+(Subscript[k, j+1 ])f((Subscript[u, j+1])(t),(Subscript[u, j])(t))-κ (Subscript[u, j])(t) . Here (Subscript[u, j])(t) is the elongation of the jth vertex, (Subscript[k, j]) is the stiffness of the jth spring connecting the vertices j-1 and j , f((Subscript[u, j-1])(t),(Subscript[u, j])(t)) is the force exerted by the j+1th vertex on the jth vertex, and κ accounts for the harmonic bounds at the original positions of the vertices along the Gosper curve. Here the forces defined.

F[k_] := force[{k, k + 1}] + force[{k, k - 1}] + force[{k}]

force[{i_, j_}] := 
Module[{p0i = vertexCoordinates[i], 
        p0j = vertexCoordinates[j], pi, pj},
        pi = p0i + {x[i], y[i]}; pj = p0j + {x[j], y[j]};
        𝓀[Max[{i, j}]] #/Sqrt[#.#]&[pj - pi]*
         (Sqrt[#.#]&[pj - pi] - Sqrt[#.#]&[p0j - p0i])]
         
(* harmonic bound to the original vertices *)
force[{k_}] := -κ {x[k], y[k]};

For small elongations, we keep only terms linear in the (Subscript[(OverHat[u]), j])={(Subscript[x, j]),(Subscript[y, j])}.

FLinearized[k_] := Normal[
   Series[F[k] /. xy:(_x | _y) :> ε xy, {ε, 0, 1}]]/ε /.
  (* endpoints *) {x[0 | 𝒹 - 1] -> 0, y[0 | 𝒹 - 1] -> 0}

This is the linearized force exerted on the second vertex.

FLinearized[2] // Collect[#, {_x, _y}]&

Forming the coefficient matrix of the {(Subscript[x, j]),(Subscript[y, j])}. In this matrix of dimension 2d×2d we first take the coefficients of the (Subscript[x, j]) and then the coefficients of the (Subscript[y, j]).

coefficientRules[𝓀_, κ_] = 
Module[{υ, xC, yC},
 (* extract x[j], y[j] and coefficient *)
 υ[p_Plus] := υ /@ (List @@ p); υ[l_List] := υ /@ l;
 υ[f_. x[j_]] := {j, f}; υ[f_. y[j_]] := {dim + j, f}; υ[0] := {};
 Flatten[Table[{xC, yC} = Collect[FLinearized[k], {_x, _y}];
               (* matrix entries *)
               {({k, #1} -> #2)& @@@ υ[xC],
               (* y[] after all the x[] *)
               ({dim + k, #1} -> #2)& @@@ υ[yC]}, {k, dim}]]];

(* form matrix from rules of nonvanishing elements *)
makeMatrix[𝓇:rules_] := 
Module[{𝓂 = Table[0, {#}, {#}]&[Max[First /@ rules]]},
      (* insert nonvanishing matrix elements *)
       Do[𝓂[[𝓇[[i, 1, 1]], 𝓇[[i, 1, 2]]]] = 
          𝓂[[𝓇[[i, 1, 1]], 𝓇[[i, 1, 2]]]] + 𝓇[[i, 2]],
          {i, Length[𝓇]}]; 𝓂]

Here are two coefficient matrices. The first M[1] uses constant spring stiffnesses and the second M[2] uses random ones with the same mean as the values of the first one (we could also use some negative spring constants [✶]) .

M[1] = makeMatrix[coefficientRules[2&, 1]];

SeedRandom[1];
M[2] = makeMatrix[
        coefficientRules[(SeedRandom[#]; 4 Random[])&, 
                         2 Random[]]];

Do[{evals[j], evecs[j]} = Eigensystem[N[M[j]]], {j, 2}]

Here are the two spectra and eigenvectors visualized. The first half are the nontrivial global modes and the second half are the superpositions of local modes due to the harmonic bounds at the lattice points and perpendicular movements. (The two intermediate modes in M[1] arise from the two endpoints.) In the case of random spring constants, we clearly see localized modes.

Do[Show[GraphicsArray[
   Block[{$DisplayFunction = Identity},
         {ListPlot[-evals[j]], ListDensityPlot[Re[evecs[j]], 
           Mesh -> False, FrameTicks -> None]}]]], {j, 2}]

The localization properties of the modes can be demonstrated by considering the overlap with the situation of all vertices are initially elongated in direction {cos(φ),sin(φ)}. The following graphics show the value of the scalar product of the modes with the initial elongation as a function of φ.

Show[GraphicsArray[
Graphics3D[{PointSize[0.002], 
Table[(* all initial positions are equal *)
      initialVector = Join[#, #]&[Table[Cos[φ], {dim}]];
      (* show overlaps *)
      MapIndexed[{Hue[#2[[1]]/dim], Point[{#2[[1]], φ, Re @ #1}]}&,
           initialVector.#& /@ evecs[#]], {φ, 0., 2.Pi, 2Pi/10}]},
     BoxRatios -> {2, 1, 1}, Axes -> True]& /@ {1, 2}]]

Now, we consider the time-dependence of an inner vertex initially elongated. The scalar products of the eigenvectors with this initial state are shown in the following two graphics.

(* one vertex is initially elongated *)
initialVector = MapAt[-1&, Table[0, {2dim}], m - 1] // N;
(* form scalar products *)
Do[coeffs[j] = (initialVector.#& /@ evecs[j]), {j, 2}];
(* show scalar products *)
Show[GraphicsArray[
 ListPlot[Re[coeffs[#]], PlotRange -> All, DisplayFunction -> Identity,
          PlotStyle -> {PointSize[0.01]}]& /@ {1, 2}]]

The time dependent elongation of all vertices follows from the superposition (Subscript[u, j])(t)=(∑_(k=1)^d)(Subscript[c, k]^((j)))(e^(i (Subscript[ω, k]) t)) where the (Subscript[c, k]^((j))) is the scalar product of the kth eigenvector with the initial elongation vector.

evI[j_][t_] := Plus @@ (Times @@@ Transpose[
               {coeffs[j] Cos[Sqrt[Abs[evals[j]]] t], evecs[j]}])

Here are the resulting time-dependent elongations shown. In the equal stiffness case, we see a constant speed propagation of the perturbation. After the two end vertices of the chain are reached, reflections occur. In the random stiffness case, the perturbation stays largely localized over the same time.

Show[GraphicsArray[Table[
(* elongation values *)
T = Table[Sqrt[Re[#.#]]& /@ Transpose[Partition[evI[j][t], dim]],
          {t, 0, 600, 2}];
(* homogeneously spaced elongation values *)
cls = #[[Round[Length[#]/2]]]& /@
         Partition[Sort[Flatten[Chop[T, 10^-12]]], dim];   
(* color function *)
cf[j][x_] = Which @@ Flatten[Append[
             MapIndexed[{x < #1, GrayLevel[(#2[[1]] - 1)/Length[cls]]}&, 
                        cls], {True, GrayLevel[1]}]];   
(* generate homogeneously colored density plot *)
ListDensityPlot[T,  PlotRange -> All, ColorFunction -> cf[j], 
                Mesh -> False, DisplayFunction -> Identity,
                MeshRange -> {{1, dim}, {0, 600}}], {j, 2}]]]

We end with an animation of the spread of the initial elongation.

frame[vec_, opts___] :=
Module[{mima, vecS, segments, elongs, colors}, 
 (* scaled elongation vectors *)
 vecS = Join[{{0, 0}}, Transpose[Partition[vec, dim]], {{0, 0}}];
 (* line along elongated vertices *)
 segments = Line /@ Partition[gosperLine + 1/2 vecS, 2, 1];
 elongs = #/Max[#]&[Sqrt[#.#]& /@ Re[vecS]]; 
 (* color segments according to elongation *)
 colors = Hue /@ (Plus @@@ Partition[elongs, 2, 1]/2);                
 (* make graphics *)
 Graphics[{{PointSize[0.03], Point[gosperLine[[m]]]},
           Transpose[{colors, segments}]}, 
           (* use uniform plot range *) PlotRange -> 
           ((# + {-1, 1})& /@ ({Min[#], Max[#]}& /@ 
            Transpose[gosperLine])), opts, 
            AspectRatio -> Automatic]];

Do[Show[GraphicsArray[Table[frame[evI[j][t]], {j, 2}]]],
   {t, 0, 600, 120}]

Do[Show[GraphicsArray[Table[frame[evI[j][t]], {j, 2}]]],
   {t, 0, 600, 5}]