Appendix

Following is some Matlab code which may be used to compute spectra of the almost Mathieu
equation, as per the paper "Reflections on the almost Mathieu operator" by Michael P.
Lamoureux. (To appear, Integral Equations and Operator Theory.)

The key peices of code are the functions Hper and Hanti, each of which returns
two tridiagonal matrices whose eigenvalues form the endpoints of the line
spectra for the almost Mathieu operator at a given rotation angle theta = p/q.
The function PlotOne plots one set of line spectra, while PlotAll puts together
a whole range of them into a pretty fractal shape.

Michael P. Lamoureux
University of Calgary


function   PlotAll()
%  Plots a range of spectra for the almost Mathieu operators, 
%  theta = p/q ranging between 0 and 1, lambda = 2.
for q = 1:25
   for p = 1:q
     if gcd(p,q) == 1
       PlotOne(p,q,2)
     end
   end
end


function   PlotOne(p,q,lambda)
%  Plots  the line spectra for one Mathieu operator at theta = p/q
%  Variables Xper, Xanti are the sorted eigenvalues of Hper, Hanti.
% In the odd q case, takes a shortcut.
[H1 H2] =  Hper(p,q,lambda); Xper  = sort([eig(H1);eig(H2)]);
if rem(q,2) = 1 
   Xanti = -Xper(q:-1:1); 
else
   [H1 H2] = Hanti(p,q,lambda); Xanti = sort([eig(H1);eig(H2)]);
end
line([Xper';Xanti'], [p/q p/q],'Color','w')


function [H1,H2] = Hper(p,q,lambda)
%      Computes two tridiagonal summands for the tridiagonal 
%      form of the Mathieu operator, periodic case.
if q == 1                     
     H1 = [2 + lambda]; H2 = [];
elseif q == 2                
     H1 = [-lambda 2; 2 lambda]; H2 = [];
elseif rem(q,2) == 0             
   q2 = q/2;
   np = rem(0:p:p*q2,q);								
   D  = lambda*cos(2*pi*np/q);        
   C  = [sqrt(2) ones(1,q2-2) sqrt(2)];  
   H1 =spdiags([ [C 0]' D' [0 C]' ], -1:1, q2+1,q2+1);
   H2 = H1(2:q2,2:q2);  
else                             
     q2 = (q+1)/2;
     np = rem(0:p:p*(q2-1),q);
     D  = lambda*cos(2*pi*np/q);        
     C  = [sqrt(2) ones(1,q2-2)];  
     H1 = spdiags([ [C 0]' D' [0 C]' ], -1:1, q2,q2);
     H2 = H1(2:q2,2:q2);
     H1(q2,q2) = H1(q2,q2) + 1;
     H2(q2-1,q2-1) = H2(q2-1,q2-1) -1;
end


function [H1,H2] = Hanti(p,q,lambda)
%      Computes two tridiagonal summands for the tridiagonal 
%      form of the Mathieu operator, antiperiodic case.
if q == 1                    
     H1 = [-lambda-2]; H2 = [];
elseif q == 2                 
     H1 = [0]; H2 = [0];
elseif rem(q,2) == 0             
     q2 = q/2;
     np = rem(p:2*p:p*(q-1),2*q);
     D  = lambda*cos(pi*np/q);        
     C  = ones(1,q2);                
     H1 = spdiags([ C' D' C' ], -1:1, q2,q2);
     H2 = H1;
     H1(1,1) = H1(1,1) + 1;  H1(q2,q2) = H1(q2,q2) - 1;
     H2(1,1) = H2(1,1) - 1;  H2(q2,q2) = H2(q2,q2) + 1;
else                             
     q2 = (q+1)/2;
     np = rem(0:p:p*(q2-1),q);
     D  = -lambda*cos(2*pi*np/q);       
     C  = [-sqrt(2) ones(1,q2-2)];  
     H1 = spdiags([ [C 0]' D' [0 C]' ], -1:1, q2,q2);
     H2 = H1(2:q2,2:q2);
     H1(q2,q2) = H1(q2,q2) - 1;
     H2(q2-1,q2-1) = H2(q2-1,q2-1) + 1;
end
    
function c = gcd(a,b)
%      Greatest Common Divisor of two integers
if  (min(a,b) <= 0)
      c = 0;
else
      c = b; r = rem(a,c);
      while (r $\sim$= 0)
             a = c; c = r; r = rem(a,c);
      end
end