%This approximates the Laplacian in 2D for functions with Neumann boundaries.
%The boundaries must be Neumann on the interval that the user chooses as the 
%   input. For example cos(x) is Neumann on [0,2*pi] while cos(pi*x) is 
%   Neumann on [0,1].
%
%This function uses the 2D discrete cosine transform. This function is in the 
%   signals toolbox of Matlab or the signal package AND image package of Octave.
%
%Sample input:
%   SpectralLaplacianNeumann2Dimensional(f,[0,1],30)
%
%The following two functions when used in the input correspond to an exact
%   solution to the laplacian that is put in the function:
%   f=@(x,y) cos(2*pi*x).*cos(2*pi*y);
%       This is the first test version. Uncomment the first For Testing area in
%       the program to see the comparison.
%   f=@(x,y) exp(cos(2*pi*x).*cos(2*pi*y));
%       This is the second test version. Uncomment the second For Testing area 
%       in the program to see the comparison.

function Result=SpectralLaplacianNeumann2Dimensional(f,Interval,n)

%%% Setting up the grid %%%

Y=0:n-1;
Y=(2*Y+1)/(2*n);
Y=Y+Interval(1);
Y=Y.*(Interval(2)-Interval(1));
Y=repmat(Y,n,1);
T=Y.';

%%% Finding the value of the function at each grid point %%%

X=f(T,Y);

%%% Finding the Laplacian %%%

X=dct2(X);  %The 2D discrete cosine transform of the evaluated function

a=repmat(0:n-1,n,1);  %A setup for the Laplacian of the interpolant functions
D=-(pi^2)*(a.^2+(a.').^2); %The coefficients of the Laplacian of the interpolant

D=D.*(1/(Interval(2)-Interval(1))^2);

%The inverse discrete cosine transform multiplied by the coefficients of the
%   interpolant's Laplacian approximates the Laplacian 
X=idct2(D.*X);

%%% Plotting the solution %%%.

figure(1);
surf(T,Y,X);
title('Approximation of 2D Laplacian','FontSize',15)

%%%%%%%%%%%For testing%%%%%%%%%%
%Uncomment the following section to compare the Laplacian of the first equation
%   given above to its approximation.

%Lf=@(x,y) (-4*pi^2*cos(2*pi*x)).*(cos(2*pi*y))+(-4*pi^2*cos(2*pi*y)).*(cos(2*pi*x));
%figure(2);
%surf(T,Y,Lf(T,Y))
%title('Actual f=@(x,y) cos(2*pi*x).*cos(2*pi*y) 2D Laplacian','FontSize',15)
%Result=max(max(abs(Lf(T,Y)-X)));

%%%%End of first testing section%%%%%%%%%%

%%%%%%%%%%For testing %%%%%%%%%%%%
%Uncomment the following section to compare the Laplacian of the second equation
%   given above to its approximation.

%Lf=@(x,y) 4*pi^2.*exp(cos(2*pi*x).*cos(2*pi*y)).*(-2*cos(2*pi*x).*cos(2*pi*y)+(sin(2*pi*y)).^2.*(cos(2*pi*y)).^2+(sin(2*pi*x)).^2.*(cos(2*pi*x)).^2);
%figure(2);
%surf(T,Y,Lf(T,Y))
%title('Actual f(x,y)= exp(cos(2*pi*x).*cos(2*pi*y)) 2D Laplacian','FontSize',15)
%Result=max(max(abs(Lf(T,Y)-X)));

%%%%%%%%%%%%% End of second testing section %%%%%%%%%%%%.