%% Second Order Finite Difference for Solving the 2D Poisson Equation on 
% a rectangle [0,1]\times[0,2] with homogeneous Neumann B.C..

m=7;
err=ones(m,1);
%% Test the accuracy
for n=1:m
dx=1/(10*n);
x=[0.5*dx:dx:1-dx/2]';
dy=dx;
y=[0.5*dy:dy:2-dy/2]';

Nx=prod(size(x)); % number of grid points in the x-direction;
Ny=prod(size(y)); % number of grid points in the y-direction;


f=-2*pi^2*cos(pi*y)*cos(pi*x'); 
uexact=cos(pi*y)*cos(pi*x')+10;    

% eigenvalues and eigenvectors for the B matrix
lambda_x=-(2*ones(Nx,1)-2*cos([0:Nx-1]'*pi*dx));
Sx=cos(x*pi*[0:Nx-1]);
iSx=inv(Sx);
lambda_y=-(2*ones(Ny,1)-2*cos([0:Ny-1]'*pi/2*dy));
Sy=cos(y*pi/2*[0:Ny-1]);
iSy=inv(Sy);
lambda=lambda_y/dy/dy*ones(1,Nx)+ones(Ny,1)*lambda_x'/dx/dx;
%S=kron(Sx,Sy);
%iS=kron(iSx, iSy);

% The eigenvector method to "invert" the B2D matrix
sol=(iSy*f*iSx')./lambda;
% notice that lambda(1,1)=0
% we choose to set 1/lambda(1,1)=0
% which will return a solution with zero sum
sol(1,1)=0; 
sol= Sy*(sol)*Sx';

% The exact solution is not unique. Since our numerical
% solution has zero sum, we should compare it to an
% exact solution with zero sum. Thus we need to subtract
% the exact solution by its sum over all indicies. 
 uexact=uexact-ones(size(uexact))*sum(sum(uexact))/Nx/Ny;

err(n)=max(max(abs(sol-uexact)));
end


close all;
hFig=figure;
set(hFig, 'Position', [0 0 500 500])
loglog(10*[1:m],err,'o'); hold all
second_order_slop=(1./(2*[1:m])').^2;
loglog(10*[1:m],second_order_slop)
xlabel('Nx');ylabel('Error')
legend('Error in max norm','Second Order Slope')

hFig=figure;
set(hFig, 'Position', [400 800 1500 400])
set(0,'DefaultTextFontSize',12,'DefaultAxesFontSize', 12)
subplot(1,4,1)
contour(x,y,f,30);colorbar;axis equal
xlabel('X');ylabel('Y')
title('The contour of function f(x,y)')
 
subplot(1,4,2)
contour(x,y,uexact,30);colorbar;axis equal
xlabel('X');ylabel('Y')
title('The contour of the exact solution u(x,y)')

subplot(1,4,3)
contour(x,y,sol,30);colorbar;axis equal
xlabel('X');ylabel('Y')
title('The contour of the numerical solution u(x,y)')

subplot(1,4,4)
contour(x,y,abs(uexact-sol),30); axis equal;colorbar
xlabel('X');ylabel('Y')
title('The contour of the error')

hFig=figure;
set(hFig, 'Position', [500 0 500 500])
surf(x,y,sol,'EdgeColor','none'); axis equal
title('The numerical solution')
xlabel('X');ylabel('Y')


