%% Second Order Finite Difference for Solving the 1D Poisson's Equation on 
%$ [0,1] with homogeneous Dirichlet B.C..
close all
clear
% Test the accuracy on m different grids
m=10;
error_2norm=ones(m,1);
error_max=ones(m,1);
% set up an exact solution 
% so the equation to be solved is
% -u''(x)= -pi*pi*sin(pi*x) with exact solution u(x)=sin(pi*x)
exact=@(x) sin(pi*x); 
exact_d1=@(x) pi*cos(pi*x);
exact_d2=@(x) -pi*pi*sin(pi*x);
for k=1:m
    n=10*k;
    % create a sparse tridiagonal matrix K
    ex = ones(n,1);
    K=spdiags([-ex 2*ex -ex], [-1 0 1], n, n);

    % create the proper uniform grid
    h=1/(n+1);
    x=[0:h:1]';x=x(2:end-1);
    % the right hand side function f(x) 
    % sampled at the grid points
    f=-exact_d2(x);
    
    % the eigenvalues of K
    lambda=2*ones(n,1)-2*cos([1:n]'*pi*h);
    % the eigenvectors of K 
    S=sin(x*pi*[1:n]);
    %% 1) A straightforward implementation of eigenvector method by
    %%    multiplying S matrix directly and its complexity is at least O(N^2)
    u= h^2* S*((S \ f)./lambda);

%     %% 2) A more efficient way (for a very large problem) is to 
%     %%    use the fact that multiplying S is equivalent to Discrete Sine Transform (O(Nlog(N)))
%     u= dst(idst(f)./lambda)*h*h;
    
    
    
    u=K\f*h*h;
    uexact=exact(x);
    error_max(k)=max(abs(u-uexact));
    error_2norm(k)=sqrt(h*sum((u-uexact).^2));
    
end

% plot the numerical and exact solutions at the finest grid

hFig=figure;
set(hFig, 'Position', [0 800 1500 400])
set(0,'DefaultTextFontSize',18,'DefaultAxesFontSize', 18)
subplot(1,2,1);
plot(x, uexact,'+'); hold all
plot(x, u,'o');
legend('exact solution','numerical solution')
axis([0 1 -0.2 1.3])
% plot the error curve and compare it with the second order slope
subplot(1,2,2);
first_order_slope=(1./(2*[1:m])').^1;
second_order_slope=(1./(2*[1:m])').^2;
loglog(2*[1:m],error_max,'+'); hold all
loglog(2*[1:m],error_2norm,'o'); hold all
loglog(2*[1:m],second_order_slope)
legend('Error in max norm','Error in  2-norm','Second Order Slope')
xlabel('Number of grid points')
ylabel('Error')