% sing_diff_mat_01.m  - analyzing singularity of the
% differentiation matrix


% Setting up the problem for a specific value of N;
  N = 2.^5;
  N = 8;
  h = 2*pi/N; x = -pi + (1:N)'*h;
  u = exp(sin(x)); uprime = cos(x).*u;

% Construct sparse second-order differentiation matrix;
  e = ones(N,1);
  D = sparse(1:N,[2:N 1],e/2.0,N,N);
  D = (D-D')/h;
    
  D
  full(D)

  
% Eigenvalues and eigenvectors of the original differentiation matrix;
  [V, E] = eig(full(D));
  [e, i] = min(abs(diag(E)));
  diag(E)
  V1 = V(:,i);
  V2 = V(:,i+1);
  
  figure(1); plot( (1:N), V1, (1:N), V2 );

  null(full(D))
  size(null(full(D)),2)
  
%  return;
  
% Removing the singularity - assume that f_1 and f_2 have to be specified
  D1 = D;
  D1(1:2,:) = 0.0; 
  D1(1,1)   = 1.0;
  D1(2,2)   = 1.0;

% The new condition number, eigenvalues and eigenvectors;  
  condest(D1)  
  [V, E] = eig(full(D1));
  diag(E)
  
  
% Dependence of the conditioning of the desingularized matrix on N;
   Nvec = 2.^(3:10);
   for N = Nvec
    h = 2*pi/N; x = -pi + (1:N)'*h;
    u = exp(sin(x)); uprime = cos(x).*u;

    % Construct sparse second-order differentiation matrix:
    e = ones(N,1);
    D =   sparse(1:N,[2:N 1],e/2.0,N,N);
    D = (D-D')/h;
    D1 = D;
    D1(1:2,:) = 0.0;
    D1(1,1)   = 1.0;
    D1(2,2)   = 1.0;
	    
    K = condest(D1);
    LL = abs(svd(full(D1)));
    S1 = min(LL);
    S2 = max(LL);
    figure(2);
    loglog(N,K,'x',N,S1,'o',N,S2,'s'), hold on
  end
  hold off
  legend('cond number','min sval','max sval')


