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

% MA/CS 615 Spring 2022, Purdue University

m=7;
err=ones(m,1);
%% Test the accuracy
for n=1:m
    
Lx=3;
Ly=2;
% the grid for K matrix    
dx=Lx/(10*n);
x=[0:dx:Lx]';
x=x(2:end-1);

% the grid for B2 matrix
dy=Ly/(10*n);
y=[0:dy:Ly]';
y=y(2:end-1);

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

% -u_xx-u_yy=f
f= (1/Lx/Lx+1/Ly/Ly)*pi^2*sin(pi*y/Ly)*sin(pi/Lx*x'); 
uexact=sin(pi*y/Ly)*sin(pi/Lx*x');    
   

lambda_x=2*ones(Nx,1)-2*cos([1:Nx]'*pi*dx/Lx);
Sx=sin(x/Lx*pi*[1:Nx]);
iSx=inv(Sx);

lambda_y=2*ones(Ny,1)-2*cos([1:Ny]'*pi*dy/Ly);
Sy=sin(y/Ly*pi*[1:Ny]);
iSy=inv(Sy);

lambda=lambda_y/dy/dy*ones(1,Nx)+ones(Ny,1)*lambda_x'/dx/dx;
 
% The eigenvector method to "invert" the matrix

f_hat=f;
 

sol=(iSy*f_hat*iSx')./lambda;
sol= Sy*(sol)*Sx';


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


close all;
hFig=figure;
set(hFig, 'Position', [0 0 500 500])
loglog(10*[1:m],err,'-.'); 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')

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