%% BVP_FDA_SimpleExample

% Finite-difference solution of 2-point ODE BVP 
%
%   d2u/dx2 = f(x) = sin(pi*x),   u(0) = 1, u(1) = 2
%
% on a grid x(j) = j/(m+1), j = 0,1,...,m,m+1.  
% using sparse tridiagonal matrix solve

%% User-specified parameters
  
  m = 9;  % Number of interior gridpoints
  alpha= 1;  % u(0)
  beta = 2;  % u(1)

  % Specify the interior grid and RHS

  h = 1/(m+1);     % grid spacing
  x = (1:m)'*h;    % mx1 column vector of interior grid points
  F = sin(pi*x);  % Specify F as a mx1 column vector
  
  F(1) = F(1) - alpha/h^2;
  F(m) = F(m) - beta/h^2;
  
  % Specify FDA 2nd derivative matrix A [m x m] on interior grid
  
  e = ones(m,1);  % m x 1 column vector of 1s
  A = spdiags([e -2*e e]/h^2, -1:1, m, m); % Specify A as sparse tridiag mx
  U = A\F; % Since A is sparse, Matlab knows to use sparse solver

  % Compute error statistics
  
  uex = -sin(pi*x)/pi^2 + alpha*(1-x) + beta*x;
  tau = A*uex - F; % Truncation error vector [m x 1]
  err = U - uex;  % Error vector
  maxerr = norm(err,inf)  % Max (infinity) norm of error
  gridnorm2_err = norm(err,2)*sqrt(h) % Grid 2-norm of error
    
%% Plot FDA vs. exact solution
  
  subplot(2,2,1)
  
  plot(x,U,'rx',[0;x;1],[alpha;uex;beta],'b-')
  xlabel('x')
  legend('FDA','Exact','Location','NorthWest')
  title(['h = ' num2str(h)])

%% Explore inverse of A and its [m+2 x m+2] augmentation Aa.  
% We index the rows/cols as i, j = 0,...,m+1 where:
%    0   corresponds to the left endpoint  (with BC U_0 = alpha) 
%    m+1 corresponds to the right endpoint (with BC U_(m+1) = beta) 

  Ainv = inv(A); 

  xa = [0; x; 1]; % Augmented col vector of gridpoints (includes endpoints)
  Fa = [alpha; F; beta]; % Augmented RHS

  % Define subdiagonal, diagonal and superdiagonal of augmented mx Aa
  % Fill in the spots we don't need to define with NaN's to show they don't 
  % matter to defining A
 
  diagAam1 = [e/h^2; 0; NaN];
  diagAa0 = [1; -2*e/h^2; 1];
  diagAap1 = [NaN; 0; e/h^2];
  Aa = spdiags([diagAam1 diagAa0 diagAap1], -1:1, m+2, m+2); % Sparse Aa 
  Aainv = inv(Aa);

  % Col j of Aainv (j=1,...,m) corresponds to h*Green's function(xa,xj)
  % Col 0 of Aainv corresponds to homogeneous solution with u(0)=1,u(1)=0
  % Col m+1 of Aainv corresponds to homogeneous solution with u(0)=0,u(1)=1
  % Also, Ainv is Aainv with the first/last row and column removed.  
  % Thus col j of Ainv (j=1,...,m) corresponds to h*Green's function(x,xj)
  % We illustrate these points in the plots below

  xj = 0.3; % Pick an interior location on the grid (otherwise arbitrary)
  j = round(xj/h);  % Index j of orresponding interior gridpoint.

% Calculate exact Green's function for xbar = xj

  Gexj= x*(xj-1); % Correct form for x <= xj
  Gexj(x>xj) = xj*(x(x>xj)-1); % Modify to correct form for x>xj
  Gexja = [0; Gexj; 0];  % Augment to include endpoint BCs.

  % Exact left and right homogeneous solutions

  Gex0a = 1-xa;
  Gexmp1a = xa;

  % Expect overall augmented finite difference solution to look like
  %
  % Ua = Gex0a*alpha + sum(j=1:m) h*Gexja*f(xj)  + Gexmp1a*beta
  %    = Gex0a*Fa(0) + sum(j=1:m) h*Gexja*F(j)  + Gexmp1a*Fa(m+1) 
  % so since Ua = Aainv*Fa, a columnwise decomposition implies
  %
  % Aainv ~= [Gex0a  h*Gexj(j=1,..,m) Gexmp1a]

  % To test this, first overplot jth col of Ainv and Aainv on h*Green's fn

  subplot(2,2,2)
  plot(xa,h*Gexja,'b-',x,Ainv(:,j),'gx',xa,Aainv(:,j+1),'ro')
  xlabel('x')
  title(['xj = ' num2str(xj) ', j = ' int2str(j)])
  legend('h*Gex(x,xj)','Ainv(:,j)','Aainv(:,j)','Location','SouthEast')

%% Lastly, interpret left and right boundary contributions to inv(Aa)

  % Overplot 0th col of Aainv on the left BC solution Gex0a(xa);

  subplot(2,2,3)
  plot(xa,Gex0a,'b-',xa,Aainv(:,1),'ro')
  xlabel('x')
  title('j=0')
  legend('Left BC uex(x)','Aainv(:,j=0)','Location','NorthEast')

  % Lastly, overplot m+1 col of Aainv on the right BC solution Gexmp1a(xa);

  subplot(2,2,4)
  plot(xa,Gexmp1a,'b-',xa,Aainv(:,m+2),'ro')
  xlabel('x')
  title('j=m+1')
  legend('Right BC uex(x)','Aainv(:,j=m+1)','Location','NorthWest')

  % For all three cases, plots show cols of Aainv (and Ainv for interior
  % points) precisely match the expectation from the exact Green's function 