%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to Problem 2 in the homework assignment,
%  Eigenvalues: Valuable Principles           
%
%  We solve the eigenvalue problem - u_{xx} - u_{yy} = lambda u
%  for (x,y) in (-1,1) x (-1,1), with u = 0 on the boundary of
%  the square.
%
%  Use the finite element method with varying meshes,
%  refining the mesh by replacing every triangle by 4.
%
%  Confirm that the convergence of eigenvalues 1, 6, 11, 16, and 21
%  is O(h^2), where h is the meshsize.
%  problem2.m       Dianne O'Leary  04/2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Modified 12-2008 to remove extra "end" statement.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  This loop creates 4 finite element meshes
%  and computes the eigenvalue approximations for each,
%  saving them in a column of lsave.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jj = 4;

for i=1:jj,
  if (i==1)
     [p,e,t]=initmesh('squareg');
  else
     [p,e,t]=refinemesh('squareg',p,e,t);
  end
  [v,l]=pdeeig('squareb1',p,e,t,1,0,1,[-inf 105]);
  sl = sort(l);
  lsave(:,i) = sl(1:21);
  npp(i) = size(p,2);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that if the iteration in pdeeig converges too slowly,
% it is possible (if the matrix is not too big) to use the
% QR algorithm to compute the eigenvalues.  Based on the code
% in pdeeig, the code would look something like this.
%
%c = 1;
%a = 0;
%d = 1;
%[p,e,t]=initmesh('squareg');
%  b = 'squareb1';
%  np=size(p,2);
%  % Boundary contributions
%  [Q,unused1,H,unused2]=assemb(b,p,e);
%  % Number of variables
%  N=size(Q,2)/np;
%  [K,M,unused3]=assema(p,t,c,a,zeros(N,1));
%  [K,unused1,B,unused2]=assempde(K,M,unused3,Q,unused1,H,unused2);
%  [unused,M]=assema(p,t,0,d,zeros(N,1));
%M = B'*M*B;
%
%zz = eig(full(K),full(M));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The true eigenvalues are lam(i,j), sorted
% and stored in lambda.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:20,
for j=1:20,
lam(i,j) = (i^2+j^2);
end
end

lambda = sort(reshape(lam,400,1))*pi^2/4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the error, print and plot it.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lamerr = lsave - lambda(1:21)*ones(1,4);

disp('Errors in the approximate eigenvalues')
disp('Lambda   Mesh1    Mesh 2    Mesh 3    Mesh 4')
disp(sprintf('  1    %7.2e  %7.2e  %7.2e  %7.2e',lamerr(1,:)))
disp(sprintf('  6    %7.2e  %7.2e  %7.2e  %7.2e',lamerr(6,:)))
disp(sprintf(' 11    %7.2e  %7.2e  %7.2e  %7.2e',lamerr(11,:)))
disp(sprintf(' 16    %7.2e  %7.2e  %7.2e  %7.2e',lamerr(16,:)))
disp(sprintf(' 21    %7.2e  %7.2e  %7.2e  %7.2e',lamerr(21,:)))

disp('Error ratios in the approximate eigenvalues')
disp('Lambda  Mesh 1vs2   Mesh 2vs3   Mesh 3vs4')
disp(sprintf('  1    %10.2e  %10.2e  %10.2e',lamerr(1,1:3)./lamerr(1,2:4)))
disp(sprintf('  6    %10.2e  %10.2e  %10.2e',lamerr(6,1:3)./lamerr(6,2:4)))
disp(sprintf(' 11    %10.2e  %10.2e  %10.2e',lamerr(11,1:3)./lamerr(11,2:4)))
disp(sprintf(' 16    %10.2e  %10.2e  %10.2e',lamerr(16,1:3)./lamerr(16,2:4)))
disp(sprintf(' 21    %10.2e  %10.2e  %10.2e',lamerr(21,1:3)./lamerr(21,2:4)))

hold off
figure(1)
loglog(npp,lamerr(1,:),npp,lamerr(6,:),npp,lamerr(11,:), ...
    npp,lamerr(16,:),npp,lamerr(21,:))
legend('\lambda_1','\lambda_6','\lambda_{11}','\lambda_{16}','\lambda_{21}')
xlabel('(approx) 1/h^2')
ylabel('error in eigenvalue')
title('Errors in eigenvalues as a function of 1/h^2')

print -depsc erreig.eps