%       program hd1dperti.m
%       Current version July 3, 2017
%       Past version: November 30, 2015
%       Computes Hausdorff dimension of the invariant set
%       obtained from iterations of the mappings:
%       theta_1(x) = (1/(3+2*lam))*x*(1 + lam*x*x*sqrt(x))
%       theta_2(x) = (1/(3+2*lam))*x*(1 + lam*x*x*sqrt(x)) + (2+lam)/(3+2*lam)
%        where 0 <= lam <= 1.
%  When lam =0, HD = ln(2)/ln(3) = 0.630929753571458
%
%     Uses the routines mat1dpert.m and thetabxlam.m
%
%  xlen = length of eigenvector
%  N+1 = number of mesh points;  Choose N to be divisable by 2
%  al is the lower matrix, aa is the approx matrix
%  mual is the largest eigenvalue of al,
%  muaa is the largest eigenvalue of aa
%  

% Input data:
       format long
       c = 0.0
       d = 1.0
       prompt = 'Input N ';  
       N = input(prompt)
       prompt = 'Input lambda ';  
       lam = input(prompt)
       blen = int16(2)
%
% Store the values of b in the vector b(l)
        b = zeros(blen,1);
        b(1) = int16(1);
        b(2) = int16(2);
        display(b)

% Calculate the mesh size
      h = (d-c)/double(N)
      xlen = N+1

%
% Index the mesh points in the domain by an index j
% of length xlen starting at 1.
%
        x = zeros(xlen,1);
        for j = 1: N +1
	    x(j) = c + double(j-1)*h;
        end
%

% For each of the mesh points x, and i = 1,2, compute the function
% aphib = (1.0/(3.0+2.0*l))*(1.0 + 3.5*lam*x*x*sqrt(x))
% which is the derivative of theta_i(x)
%
          aphib = zeros(xlen,2);
          for k = 1: xlen
          aphib(k,1) = (1.0/(3.0+2.0*lam))*(1.0 + 3.5*lam*x(k)*x(k)*sqrt(x(k)));
          aphib(k,2) = aphib(k,1);
          end
%
%  For each of the points thetab(x) = thetabx, 
%  find the index number ix of the left mesh point
%  of the subinterval to which it belongs.
%  Note that we must decide how to choose the subinterval 
%  when the points lie on the boundary of two subintervals
%
ix = zeros(xlen,blen);    
Lx = ones(xlen,blen,2);
for k = 1:xlen
    for l = 1:blen
       thetabxkl = thetabxlam(x(k),b(l),lam);
       ix(k,l) = floor((thetabxkl-c)/h) + 1;
%
%  If the index is outside the allowable limits, print the results
         if (ix(k,l) < 1)
	       display(['(ix(k,l) < 1', k, l, ix(k,l), thetabxkl])
         end
         if (ix(k,l) > xlen)
	       display(['ix(k,l) > xlen', k, l, ix(k,l), thetabxkl])
         end
% If the index is equal to xlen, then reset to one less
         if (ix(k,l) == xlen)
	   ix(k,l) = ix(k,l) -1;
         end
%   We next compute the Lagrange basis functions Lx corresponding 
%   to each of these degrees of freedom at the point thetabx
%
    ixkl = ix(k,l);
    Lx(k,l,2) = (thetabxkl - x(ixkl))/h;
    Lx(k,l,1) = 1.0 - Lx(k,l,2);
  end
end
%
%
% Compute the root of log(max evalue(s)) = 0 using the Matlab routine fzero.
%
% First do the computation for the approximate matrix
mtype = 0;
funa = @(s) log(eigs(mat1dpert(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx, lam),1));
sa = fzero(funa, [0.01, 0.99]);
muaa = eigs(mat1dpert(mtype, sa, h, ix, x, b, aphib, xlen, blen, Lx,lam),1);
muaa1 = muaa - 1;
disp ('Approximation of Hausdorff dimension s.')
sa
disp ('The eigenvalue of the  approximate matrix minus one.')
muaa1
%
% Next do the computation for the lower matrix
mtype = -1;
funl = @(s) log(eigs(mat1dpert(mtype, s, h, ix, x, b, aphib, xlen, blen, Lx,lam),1));
sl = fzero(funl, sa);
As = mat1dpert(mtype, sl, h, ix, x, b, aphib, xlen, blen, Lx,lam);
[val, mual] = eigs(As,1);
mual1 = mual - 1;
% Add loop to ensure mual > 0  and min(As*val - val) >=0
% First make sure that val is a strictly positive eigenvector
if (max(val) <= 0) val = -val;
end
minval = min(As*val - val);
while  (mual1 <= 0 || minval < 0)
   sl = sl - eps;
As = mat1dpert(mtype, sl, h, ix, x, b, aphib, xlen, blen, Lx,lam);
   [val, mual] = eigs(As,1);
   mual1 = mual - 1;
   if (max(val) <= 0) val = -val;
   end
minval =   min(As*val - val);
end
disp ('Lower bound on the Hausdorff dimension s.')
sl
disp ('The eigenvalue of the  lower matrix minus one.')
% disp ('The value should be >=0.')
mual1
% The following value should be >=0
minval
%
%       end program hd1dperti.m
