CHAPTER 16, Two-Dimensional Substitutions.

16.1 Substitution Sequences.

We look at certain two-dimensional substitutions and their tilings.

PlotResults[x_]:=MatrixPlot[x,ColorRules-> {0->White,1-> Black,2-> Red},Mesh->All]

θ1[x_]:=
Switch[x,
(* To change our transformation θ, do it here! *)
0,{{1,2},{0,1}},1,{{2,0},{0,1}},2,{{2,0},{1,2}},_,Message[θ::ArgumentOutOfRange,x];Abort[]
]


θ[x_]:=
If[ArrayDepth[x]==0,
(* x is number *)
θ1[x],
(* Else *)
(* x is matrix *)
ArrayFlatten[Map[θ1,x,{ArrayDepth[x]}]]
]

TPow[x_,power_]:=(
res=x;
For[i=1,i≤power,i++,
res=θ[res];
];
res
)

Example 1.

We are looking at the two-dimensional substitution:
θ (0) =  (({{1, 2}, {0, 1}})), θ (1) =  (({{2, 0}, {0, 1}})), θ (2) =  (({{2, 0}, {1, 2}})),

θ[0]//MatrixForm
θ[1]//MatrixForm
θ[2]//MatrixForm

(* Applying (θ^2)[0] *)
θ[θ[0]]//MatrixForm
TPow[0,3]//MatrixForm

PlotResults[TPow[0,3]]

Example 2. A 3×3 Substitution.

θ1[x_]:=
Switch[x,0,{{1,2,0},{0,1,2},{0,2,1}},
1,{{1,2,0},{2,0,1},{0,1,2}},2,{{1,0,2},{2,1,0},{1,2,0}},
_,Message[θ::ArgumentOutOfRange,x];Abort[]
]

θ[0]//MatrixForm
θ[1]//MatrixForm
θ[2]//MatrixForm

θ[θ[0]]//MatrixForm
TPow[0,2]//MatrixForm

PlotResults[
TPow[0,2]
]

Example 3. A 4×4 Substitution.

θ1[x_]:=
Switch[x,
0,{{1,2,0,1},{0,1,2,0},{0,2,1,2},{2,1,0,1}},
1,{{1,2,0,1},{2,0,1,2},{0,1,2,0},{1,0,1,2}},
2,{{1,0,2,0},{2,1,0,1},{1,2,0,2},{0,2,0,1}},
_,Message[θ::ArgumentOutOfRange,x];Abort[]
]

θ[0]//MatrixForm
θ[1]//MatrixForm
θ[2]//MatrixForm

θ[θ[0]]//MatrixForm
TPow[0,2]//MatrixForm

PlotResults[
TPow[0,2]
]

16.2 The Rauzy Substitution and Incidence Matrices.

Let us recall how the Fibonacci sequence F(n), n=0, 1, 2, …, is constructed recursively: F(0)=0, F(1)=1, and F(n+2)=F(n+1)+F(n), for n≥0.

F[0]=0; F[1]=1; 
Do[F[i+2]=F[i+1]+F[i],{i,0,20}]
Table[F[i],{i,0,20}]

In a similar way, the Tribonacci sequence T(n), n=0, 1, 2, …, is constructed: T(0)=0, T(1)=0, T(2)=1, and T(n+2)=T(n+1)+T(n)+T(n-1), for n≥1.

T[0]=0; T[1]=0; T[2]=1 
Do[T[i+2]=T[i+1]+T[i]+T[i-1],{i,1,20}]
Table[T[i],{i,0,20}]

The matrix of the Fibonacci substitution is given by.:

A={{1,1},{1,0}};MatrixForm[A]
Det[A-z*IdentityMatrix[2] ]
Eigenvalues[A]

Doing the same for the Tribonacci substitution:

T={{1,1,1},{1,0,0},{0,1,0}};MatrixForm[T]
Det[T-z*IdentityMatrix[3] ]

Consequently, the eigenvalues of T are given by:

Solve[1+z+z^2-z^3==0,z]

i.e., ρ=(1/3) (1+(((19-3 (Sqrt[33])))^(1/3))+(((19+3 (Sqrt[33])))^(1/3)))=1.83928… and λ=(1/3)-(1/6) (1+i (Sqrt[3])) (((19-3 (Sqrt[33])))^(1/3))-(1/6) (1-i (Sqrt[3])) (((19+3 (Sqrt[33])))^(1/3))=-0.419643+0.606291 i, (approximately), together with the conjugate of λ. Note that the matrix T is primitive (the entries of (T^n), for some n, are all positive). Recall that we use A.B for matrix multiplication. In particular, ρ is a Pisot number.

MatrixForm[T.T.T]

16.2 Exercises.

1.   Generalize the Fibonacci and Tribonacci substitutions in a natural way to find a substitution on A={1,2,3,4}, whose incidence matrix has eigenvalues satisfying the characteristic equation 1+z+(z^2)+(z^3)-(z^4)=0:

16.3 Rauzy Fractals.

Clear[f,dlst,pt,cr,ptlst]
dlst=Table[Random[Integer,{1,3}],{n,250000}];
in=N[{{1,0},{-1/2,Sqrt[3]/2},{-1/2,-Sqrt[3]/2}}];
f[j_,{x_,y_}]:={-y,x}/N[Sqrt[3]]+in[[j]]
pt={0.5,0.5};
cr[n_]=If[n-1==0,RGBColor[0,0,1],If[n-2==0,RGBColor[0,1,0],If[n-3==0,RGBColor[1,0,0],RGBColor[0,0,0]]]];
ptlst=Table[{cr[dlst[[j]]],Point[pt=f[dlst[[j]],Sequence[pt]]]},{j,Length[dlst]}];
Show[Graphics[Join[{PointSize[.001]},ptlst]],AspectRatio->Automatic,PlotRange->All]

Clear[f,dlst,pt,cr,ptlst]
dlst=Table[Random[Integer,{1,3}],{n,600000}];
rotate[theta_]:={{Cos[theta],-Sin[theta]},{Sin[theta],Cos[theta]}};
NSolve[x^3-x^2-x-1==0,x]
an=Table[Arg[x^n]/.NSolve[x^3-x^2-x-1==0,x][[1]],{n,1,3}];
M=Table[(Abs[x^n]/.NSolve[x^3-x^2-x-1==0,x][[1]])*rotate[an[[n]]],{n,1,3}];
in=Table[{Re[x^n]/.NSolve[x^3-x^2-x-1==0,x][[1]],Im[x^n]/.NSolve[x^3-x^2-x-1==0,x][[1]]},{n,1,3}];
f[j_,{x_,y_}]:=M[[j]].{x,y}+in[[j]]

pt={0.5,0.5};

cr[n_]=If[n-1==0,RGBColor[0,0,1],If[n-2==0,RGBColor[0,1,0],If[n-3==0,RGBColor[1,0,0],RGBColor[0,0,0]]]];
ptlst=Table[{cr[dlst[[j]]],Point[pt=f[dlst[[j]],Sequence[pt]]]},{j,Length[dlst]}];

Show[Graphics[Join[{PointSize[.001]},ptlst]],AspectRatio->Automatic,PlotRange->All]

16.4 More Rauzy Fractals.

Clear[f,dlst,pt,cr,ptlst]
dlst=Table[Random[Integer,{1,3}],{n,200000}];
rotate[theta_]:={{Cos[theta],-Sin[theta]},{Sin[theta],Cos[theta]}};
NSolve[x^3-x^2-1==0,x]
an=Table[Arg[x^n]/.NSolve[x^3-x^2-1==0,x][[1]],{n,1,3}];
M=Table[(Abs[x^n]/.NSolve[x^3-x^2-1==0,x][[1]])*rotate[an[[n]]],{n,1,3}];
in=Table[{Re[x^n]/.NSolve[x^3-x^2-1==0,x][[1]],Im[x^n]/.NSolve[x^3-x^2-1==0,x][[1]]},{n,1,3}];
f[j_,{x_,y_}]:=M[[j]].{x,y}+in[[j]]

pt={0.5,0.5};

cr[n_]=If[n-1==0,RGBColor[0,0,1],If[n-2==0,RGBColor[0,1,0],If[n-3==0,RGBColor[1,0,0],RGBColor[0,0,0]]]];
ptlst=Table[{cr[dlst[[j]]],Point[pt=f[dlst[[j]],Sequence[pt]]]},{j,Length[dlst]}];

Show[Graphics[Join[{PointSize[.001]},ptlst]],AspectRatio->Automatic,PlotRange->All]

Module[{O=12,iter=5,
w=(2-(1+I*Sqrt[3])(19-3Sqrt[33])^(1/3)+I(I+Sqrt[3])(19+3Sqrt[33])^(1/3))/6 //N,
step, rauzyPoints, toThreeparts},
step[{a___,b_,c_}]:=
If[b==c==1,{{a,b,c,0}},{{a,b,c,0},{a,b,c,1}}];
rauzyPoints=Function[ws,(ws.#)&/@Nest[
Flatten[step[#]&/@#,1]&,{{0,0}},O]][w^Range[0,O+1]];
toThreeParts[rp_]:={(w#)&/@rp,(w^2#+w^2)&/@rp,
(w^3#+w^3+w^2)&/@rp};
Show[Graphics[{MapIndexed[{
Hue[Sin[#2[[1]]]],PointSize[0.002],
Point[{Im[#],Re[#]}]&/@#}&,
Nest[Flatten[toThreeParts/@#,1]&,
{rauzyPoints},iter],{1}]}],
AspectRatio->Automatic]]

16.5 Dragon Curves.

Directions="1";(*Directions pattern, 1 - turn left, 2 - turn right *)
Edge=1;
Angle=Pi/4;

Iterate=Input["Please enter the level of the dragon curve"]; 
If[Iterate>15,Iterate=15];

Do[
CurrDirection=Directions; 
NewDirectionOne=StringJoin[CurrDirection,"1"];
Dirlength=(StringLength[Directions]-1)/2;
NewDirectionTwo=StringJoin[StringTake[Directions,Dirlength],"0",StringTake[Directions,-1*Dirlength]];
Directions=StringJoin[NewDirectionOne,NewDirectionTwo];
Angle=Angle+Pi/4;
If[Angle>=Pi*2,Angle=0];
Edge=Edge/Sqrt[2]; 
,{n,2,Iterate}];


X=0;
Y=0;
ListX=Table[0*n,{n,2^Iterate+1}];
ListY=Table[0*n,{n,2^Iterate+1}];
ListX[[1]]=X;
ListY[[1]]=Y;
X=X+Cos[Angle]*Edge;
Y=Y+Sin[Angle]*Edge;
ListX[[2]]=X;
ListY[[2]]=Y;


Do[
DirectionStep=StringTake[Directions,{n-2}];
If[StringMatchQ[DirectionStep,"1"],Angle=Angle-Pi/2,Angle=Angle+Pi/2];
X=X+Cos[Angle]*Edge;
Y=Y+Sin[Angle]*Edge;

ListX[[n]]=X;
ListY[[n]]=Y;
,{n,3,2^Iterate+1}];

Pairs=Table[{ListX[[n]],ListY[[n]]},{n,2^Iterate+1}];
ListPlot[Pairs,Joined->True,AspectRatio->Automatic,Axes->False]