Localization on a Harmonic Chain
with Zero Displacement at the Boundaries

Amy Kolan, St. Olaf College, Northfield, MN, 55057; kolan@stolaf.edu.  
Barry Cipra, 305 Oxford St., Northfield, MN 55057; cipra@stolaf.edu.  
Bill Titus, Carleton College, Northfield, MN 55057; btitus@physics.carleton.edu

Introduction

This Mathematica notebook contains the code used to generate the figures in the article "Exploring Localization in Nonperiodic Systems" which appears in the June 1995 issue of Computers in Physics.  Please refer to that article for a discussion of localization and the notation used in this notebook.  This notebook is not meant to be self-contained, but to be used in conjunction with the article.

Some Basic Function Definitions

Fundamental transfer matrix (See Eq. (11))

Here  m  is the mass,  k  is  the spring constant, and   w   is   ω2.

Clear[m, w, k, T]
T[m_, w_, k_] = {{1 - m*w/k, 1/k}, {-m*w, 1}};

Local eigenvalue functions

This function is used to determine  α  from the trace of a transfer matrix (see Eq. (16)).

Clear[alpha, g]
alpha[g_] := (1/2)*(g[[1, 1]] + g[[2, 2]])

This function determines the log of the magnitude of the largest local eigenvalue for the transfer matrix when  |α| >= 1  (see Eq. (17)).  For  |α| < 1, the log of the eigenvalue is imaginary and, for plotting purpose, the result is set to zero.

Clear[eta, alf]
eta[alf_ /; Abs[alf] >= 1] := ArcCosh[Abs[alf]]
eta[alf_ /; Abs[alf] < 1]  := 0

Global eigenvalues:  the characteristic equation for  w

The eigenvalues for  w  are those that cause  qN+1[w]  to be zero.  From Eq. (15),  qN+1[w]  is equal to the  {1,2}  element of the product of the transfer matrices across the chain.   In the function below,   n  corresponds to  N.

Clear[m, k, w, n, qFinal]
qFinal[w_, n_] :=
  Module[
    {temp, i}, 
    temp = T[m[0], w, k[0]] // Expand;
    For[i = 1, i <= n, i++,
        temp = T[m[i], w, k[i]] . temp // Expand];
    temp[[1,2]]
    ];

A Specific Run

The labels for the masses run from  i = 0  to  i = n+1,  where the first and last displacements are zero.  The labels for the spring constants go from  i = 0  to  i = n .  The results for  w  and  qi[w]  are independent of the value for  m[0].

Values for the parameters  m  and  k

Clear[imax, m, k]

n = 20;

For[i = 0, i <= n, i++, k[i] = 1];

For[i = 1, i <= n+1, i++, m[i] = 4];

   m[0]   = 0;
   m[7]   = 1;

Plot of  qN+1[w]   as a function of  w.

The normal modes occur when  qN+1[w] = 0.  The following plot visually shows the distribution of the zeroes of   qN+1[w] and hence the allowed values of  w.  Remember,  w = ω2.

Clear[w]
Plot[Evaluate[qFinal[w, n]], {w, 0.0, 2.5},
     Frame      -> True,
     PlotStyle  -> Thickness[0.00],
     PlotRange  -> {-5, 5},
     FrameLabel -> {"w", "qfinal"},
     PlotPoints -> 50];

Numerical solution for eigenvalues using  NSolve[]

Here   NSolve[]  is used to numerically determine the allowed values of  w.

Clear[w, equ, list1, wList]
equ = qFinal[w, n] // Expand;
list1 = NSolve[equ == 0, w, WorkingPrecision -> 64];
wList = w /. list1 // N[#, 16]&

Here's a quick check to see that the values of  w  do indeed cause  qFinal[w] = 0.  Note that the value of  qFinal  for the last eigenvalue is not zero, indicating an extreme sensitivity in the determination of the eigenvalue, i.e., for this eigenvalue, a small change in  w  will produce a large change in  qFinal.

Table[qFinal[wList[[i]], n], {i, 1, Length[wList]}]

Various Displays of Position and Local Eigenvalues (Figures 3 - 5)

wNewList = Take[wList, {n-1, n}]

Clear[qPlot, qSqrPlot, qCumPlot, logqSqrPlot, etaPlot]
For[j = 1, j <= Length[wNewList], j++,

    w = wNewList[[j]];  
    
    Clear[t];
    t[0] = T[m[0], w, k[0]];
    For[i = 1, i <= n, i++,
        t[i] = T[m[i], w, k[i]] . t[i-1]];
            
    Clear[Q];
    Q[0] = {0, 1};
    For[i = 1, i <= n+1, i++,
        Q[i] = t[i-1] . Q[0]];
                 
    Clear[qSumNorm];
    qSumNorm = 0;
    For[i = 1, i <= n, i++,
        qSumNorm = qSumNorm + m[i] Q[i][[1]]^2;
       ];  
    
    Clear[q];
    For[i = 0, i <= n+1, i++,
        q[i] = Sqrt[1/qSumNorm]*Q[i][[1]]];
             
    Clear[qCumNorm];
    qCumNorm[0] = 0;
    For[i = 1, i <= n, i++,
        qCumNorm[i] = qCumNorm[i-1] + m[i] q[i]^2];
    qCumNorm[n+1] = 1.0;
 
    temp1 =  
      ListPlot[Table[{i, q[i]}, {i, 0, n+1}],
               PlotStyle  -> {Thickness[0.00], PointSize[0.02]},
               DisplayFunction -> Identity];
                              
    temp2 =  
      ListPlot[Table[{i, q[i]}, {i, 0, n+1}],
               PlotStyle  -> Thickness[0.00],
               PlotJoined -> True,
               DisplayFunction -> Identity];
    
    qPlot[j] =
      Show[{temp1, temp2},
           Frame      -> True,
           FrameLabel -> {"j", "q", 
              FontForm[
                  StringJoin["w = ", 
                    ToString[StringForm["``", 
                             NumberForm[Sqrt[w], {5, 6},
                             NumberPadding -> {"", "0"}]]]],
                  {"Symbol", 10}], ""}, 
           AxesOrigin -> {0, 0},            
           PlotRange  -> {{-1, n+2}, All},
           FrameTicks -> {{1, 5, 10, 15, 20}, Automatic}]; 
                                   
   temp1 =  
      ListPlot[Table[{i, q[i]^2}, {i, 0, n+1}],
               PlotStyle  -> {Thickness[0.00], PointSize[0.02]},
               DisplayFunction -> Identity];
                              
    temp2 =  
      ListPlot[Table[{i, q[i]^2}, {i, 0, n+1}],
               PlotStyle  -> Thickness[0.00],
               PlotJoined -> True,
               DisplayFunction -> Identity];
     
    qSqrPlot[j] = 
      Show[{temp1, temp2},
           Frame      -> True,
           FrameLabel -> {"j", "qSqr", 
              FontForm[
                  StringJoin["w = ", 
                    ToString[StringForm["``", 
                             NumberForm[Sqrt[w], {5, 6},
                             NumberPadding -> {"", "0"}]]]],
                  {"Symbol", 10}], ""}, 
           AxesOrigin -> {0, 0},
           FrameTicks -> {{1, 5, 10, 15, 20}, Automatic},
           PlotRange  -> {{-1, n+2}, All}];
   
   temp1 =  
      ListPlot[Table[{i, qCumNorm[i]}, {i, 0, n+1}],
               PlotStyle  -> {Thickness[0.00], PointSize[0.02]},
               DisplayFunction -> Identity];
                              
    temp2 =  
      ListPlot[Table[{i, qCumNorm[i]}, {i, 0, n+1}],
               PlotStyle  -> Thickness[0.00],
               PlotJoined -> True,
               DisplayFunction -> Identity];
          
    qCumPlot[j] = 
      Show[{temp1, temp2},
           Frame      -> True,
           FrameLabel -> {"j", "qCum", 
              FontForm[
                  StringJoin["w = ", 
                    ToString[StringForm["``", 
                             NumberForm[Sqrt[w], {5, 6},
                             NumberPadding -> {"", "0"}]]]],
                  {"Symbol", 10}], ""}, 
            AxesOrigin -> {0, 0},
            FrameTicks -> {{1, 5, 10, 15, 20}, Automatic},
            PlotRange  -> {{-1, n+2}, All}];
    
    temp1 =  
      ListPlot[Table[{i, Log[Abs[q[i]]]}, {i, 1, n}],
               PlotStyle  -> {Thickness[0.00], PointSize[0.02]},
               DisplayFunction -> Identity];
                              
    temp2 =  
      ListPlot[Table[{i, Log[Abs[q[i]]]}, {i, 1, n}],
               PlotStyle  -> Thickness[0.00],
               PlotJoined -> True,
               DisplayFunction -> Identity];
           
    logqSqrPlot[j] = 
      Show[{temp1, temp2},
           Frame      -> True,
           FrameLabel -> {"j", "Log[|q|]", 
              FontForm[
                  StringJoin["w = ", 
                    ToString[StringForm["``", 
                             NumberForm[Sqrt[w], {5, 6},
                             NumberPadding -> {"", "0"}]]]],
                  {"Symbol", 10}], ""}, 
           AxesOrigin -> {0, 0},
           FrameTicks -> {{1, 5, 10, 15, 20}, Automatic},
           PlotRange  -> {{-1, n+2}, All},
           DisplayFunction -> Identity];
    
    temp1 =  
      ListPlot[Table[{i, eta[alpha[t[i]]]}, {i, 0, n}],
               PlotStyle  -> {Thickness[0.00], PointSize[0.02]},
               DisplayFunction -> Identity];
    temp2 =  
      ListPlot[Table[{i, eta[alpha[t[i]]]}, {i, 0, n}],
               PlotStyle  -> Thickness[0.00],
               PlotJoined -> True,
               DisplayFunction -> Identity];
 
     etaPlot[j] = 
       Show[{temp1, temp2},
            Frame      -> True,
            FrameLabel -> {"j", "eta", 
              FontForm[
                  StringJoin["w = ", 
                    ToString[StringForm["``", 
                             NumberForm[Sqrt[w], {5, 6},
                             NumberPadding -> {"", "0"}]]]],
                  {"Symbol", 10}], ""}, 
             AxesOrigin -> None,
             FrameTicks -> {{1, 5, 10, 15, 20}, Automatic},
             PlotRange  -> {{-1, n+2}, All}]            
  ];

Show[GraphicsArray[{qPlot[1], qPlot[2]}]];

Show[GraphicsArray[{{qSqrPlot[1], qSqrPlot[2]},
                    {qCumPlot[1], qCumPlot[2]}}]];

Show[GraphicsArray[{{logqSqrPlot[1], logqSqrPlot[2]},
                    {etaPlot[1], etaPlot[2]}}]];