% This function is to help with homework problem 5.4 in MATH 393C.
% It contains the following functions:
%
% DRIVER_singletest: A spectral collocation method for solving a 2 point
% boundary value problem. 
%
% DRIVER_box_exact_soln: Spectral collocation method for solving the
% Helmholtz equation on a box. We solve a problem with a known analytic
% solution.
%
% DRIVER_box: Spectral collocation method for solving the Laplace 
% equation on a box for general Dirichlet data.

function HW05p4

DRIVER_singletest
DRIVER_box_exact_soln
DRIVER_box

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function solves a 2 points boundary value problem using a spectral
% collocation method with a known solution.
%
% This function demonstrates that for a spectral method, it is misleading
% to look at the linear interpolant between the nodes. The nodes are sparse
% enough that the linear interpolant may look completely wrong even though
% each spot is extremely accurate.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_singletest

%%% Pick the number of discretization nodes, generate the Chebyshev grid
%%% and the spectral differentiation operator.
n      = 21;
[D,xx] = LOCAL_cheb(n-1);
xx     = xx(end:(-1):1);
D      = D(end:(-1):1,end:(-1):1);

%%% Build the elliptic operator.
kh     = 12;
A      = -D*D - kh*kh*eye(n);

%%% Partition the grid.
Iext = [1,n];
Iint = 2:(n-1);

%%% Set the exact solution:
%%% We also generate a more finely sampled solution for plotting purposes.
xxfine = linspace(-1,1,1000)';
ufine  = cos(kh*xxfine);
u      = cos(kh*xx);

%%% Solve the equation
u_approx       = zeros(n,1);
u_approx(Iext) = u(Iext); % Use the given Dirichlet data.
u_approx(Iint) = -A(Iint,Iint)\(A(Iint,Iext)*u_approx(Iext)); % Solve.
fprintf(1,'Max error          = %12.3e\n',max(abs(u - u_approx)))
fprintf(1,'cond(A(Iint,Iint)) = %12.3e\n',cond(A(Iint,Iint)))

figure(1)
plot(xxfine,ufine,'b',xx,u_approx,'r.-')
legend('Exact solution','Computed solution')

figure(2)
[u_approx_interp]=barylag([xx,u_approx],xxfine);
plot(xxfine,ufine,'b',...
     xxfine,u_approx_interp,'g--',...
     xx,u_approx,'r.','MarkerSize',15)
legend('Exact solution','Computed solution (interpolated)')

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_box_exact_soln

n  = 30;
kh = 10;
[L,D1,D2,xvec] = LOCAL_get_L_and_D(n,1);
A = -L - kh*kh*eye(n*n);

[XX1,XX2] = meshgrid(xvec);
xx = [reshape(XX1,1,numel(XX1));...
      reshape(XX2,1,numel(XX2))];
hmin = min(abs(xvec(2) - xvec(1)));

Iint = find((abs(xx(1,:)) < (1 - 0.5*hmin)) & ...
            (abs(xx(2,:)) < (1 - 0.5*hmin)));
Iext = find((abs(xx(1,:)) > (1 - 0.5*hmin)) | ...
            (abs(xx(2,:)) > (1 - 0.5*hmin)));

%%% Set an exact solution.
dist           = sqrt((xx(1,:) - 1.5).^2 + (xx(2,:) - 0.25).^2)';
u              = bessely(0,kh*dist);

%%% Compute the approximate solution.
u_approx       = zeros(size(u));
u_approx(Iext) = u(Iext);
u_approx(Iint) = -A(Iint,Iint)\(A(Iint,Iext)*u_approx(Iext));
fprintf(1,'Max error          = %12.3e\n',max(abs(u - u_approx)))
fprintf(1,'cond(A(Iint,Iint)) = %12.3e\n',cond(A(Iint,Iint)))
        
%%% Plot the solution and the error.
figure(1)
subplot(1,2,1)
mesh(XX1,XX2,reshape(u_approx,n,n))
title('Computed solution')

subplot(1,2,2)
mesh(XX1,XX2,reshape(u - u_approx,n,n))
title('Error')

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function DRIVER_box

m  = 30;
n  = 2*m+1;
[L,D1,D2,xvec] = LOCAL_get_L_and_D(n,1);
A  = -L;

[XX1,XX2] = meshgrid(xvec);
xx = [reshape(XX1,1,numel(XX1));...
      reshape(XX2,1,numel(XX2))];
hmin = min(abs(xvec(2) - xvec(1)));

Iint = find((abs(xx(1,:)) < (1 - 0.5*hmin)) & ...
            (abs(xx(2,:)) < (1 - 0.5*hmin)));
Iext = find((abs(xx(1,:)) > (1 - 0.5*hmin)) | ...
            (abs(xx(2,:)) > (1 - 0.5*hmin)));

%%% Set a boundary condition.
u_approx       = zeros(n*n,1);
%u_approx(Iext) = (cos(11*xx(2,Iext)).*exp(xx(1,Iext)))';
u_approx(Iext) = (((xx(1,Iext)-0.5).^2 + (xx(2,Iext)-1).^2).^0.25)';
        
%%% Compute the approximate solution.
u_approx(Iint) = -A(Iint,Iint)\(A(Iint,Iext)*u_approx(Iext));
%fprintf(1,'cond(A(Iint,Iint)) = %12.3e\n',cond(A(Iint,Iint)))
fprintf(1,'uapprox(0,0) = %20.15f\n',u_approx(floor(n*n/2)+1))
        
%%% Plot the solution and the error.
figure(1)
subplot(1,1,1)
mesh(XX1,XX2,reshape(u_approx,n,n))
title('Computed solution')

keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [p]=barylag(data,x)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% barylag.m
%
% Interpolates the given data using the Barycentric
% Lagrange Interpolation formula. Vectorized to remove all loops
%
% data - a two column vector where column one contains the
%        nodes and column two contains the function value 
%        at the nodes
% p - interpolated data. Column one is just the 
%     fine mesh x, and column two is interpolated data 
%
% Reference:
%
% (1) Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange 
%     Interpolation" 
%     http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/berrut.ps.gz
% (2) Walter Gaustschi, "Numerical Analysis, An Introduction" (1997) pp. 94-95
%
%
% Written by: Greg von Winckel       03/07/04
% Contact:    gregvw@chtm.unm.edu 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=length(data);     N=length(x);

% Compute the barycentric weights
X=repmat(data(:,1),1,M);

% matrix of weights
W=repmat(1./prod(X-X.'+eye(M),1),N,1);

% Get distances between nodes and interpolation points
xdist=repmat(x,1,M)-repmat(data(:,1).',N,1);

% Find all of the elements where the interpolation point is on a node
[fixi,fixj]=find(xdist==0);

% Use NaNs as a place-holder
xdist(fixi,fixj)=NaN;
H=W./xdist;

% Compute the interpolated polynomial
p=(H*data(:,2))./sum(H,2);

% Replace NaNs with the given exact values. 
p(fixi)=data(fixj,2);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function computes the Chebyshev nodes on [-1,1].
% It also computes a differentiation operator.
% It is modified from a code by Nick Trefethen.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [D,x] = LOCAL_cheb(N)

if N==0
  D=0; 
  x=1; 
  return
end
x = cos(pi*(0:N)/N)'; 
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';                  
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D'));                 % diagonal entries

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function constructs the spectral differentiation matrices.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [L,D1,D2,xvec] = LOCAL_get_L_and_D(N_side,a)

[D,xvec] = LOCAL_cheb(N_side-1);
xvec     = a*xvec(end:(-1):1);
D        = (1/a)*D;
I        = eye(N_side);
D1       = -kron(D,I);
D2       = -kron(I,D);
Dsq      = D^2; 
L        = kron(I,Dsq) + kron(Dsq,I);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
