% lap5pt.m  -- solve Dirichlet's problem for the Poisson equation
%LAP5PT  5 point Laplacian difference method for Poisson equation
%
%         Laplacian_h u = f  on interior grid points
%                     u = g  on boundary grid points
%
% using 5-point Laplacian.  Domain is unit square.
% Douglas N. Arnold, adapted from poisson.m of Randy Leveque

% supply f and, if known, exact solution u as functions
%% example 0: u a polynomial of < degree 4 in x and y (so method should be exact at nodes)
%u = @(x,y) x.^3.*(y.^3 - 3.*y.^2 + 2)/2 - 100*(x.^2-x).*y;
%f = @(x,y) 3*x.*y.^3 - 9*x.*y.^2 + 3*x.^3.*y -3*x.^3 - 200*y + 6*x;
%% example 1: u = sin pi x  sin pi y
%in(pi*x).*sin(pi*y);
%f = @(x,y) -2*pi^2*sin(pi*x).*sin(pi*y);
%% example 2:
%u = @(x,y) exp(-50.*(x-2/3).^2-200.*(y-1/3).^2)/2 - cos(9/2*pi*x)/8;
%f = @(x,y) exp(-50.*(x-2/3).^2-200.*(y-1/3).^2)*250.*(391-240*x+180*x.^2-1920*y+2880*y.^2)/9 +81*pi^2*cos(9/2*pi*x)/32;
%% example 3:
u = @(x,y) .02*sin(8*pi*x*y)  % boundary data
f = @(x,y) zeros(size(x)) % solve homogeneous Laplace equation

% set mesh size
n = 128;
h = 1/n;
m = n-1;
x = linspace(0,1,n+1);   % x grid points including boundaries
y = linspace(0,1,n+1);   % y grid points including boundaries

[X,Y] = meshgrid(x,y);      % 2d arrays of x,y values
X = X';                     % transpose so that X(i,j),Y(i,j) are
Y = Y';                     % coordinates of (i,j) grid point, i.e., ((i-1)h,(j-1)h)

Iint = 2:n;                 % indices of interior points in x
Jint = 2:n;                 % indices of interior points in y
Xint = X(Iint,Jint);        % interior points
Yint = Y(Iint,Jint);

fh = f(Xint,Yint);          % fh is a mesh function on interior grid points

% set uh on the boundary to given values (interior values will be ignored)
uh = u(X,Y);        % can use true solution if it is known

% form 5-point Laplacian matrix (times h^2)
I = speye(m);
e = ones(m,1);
A = spdiags([e,-4*e,e],[-1,0,1],m,m);
J = spdiags([e,e],[-1,1],m,m);
Lh = kron(I,A) + kron(J,I);

% multiply right hand side by h^2 and adjust for boundary condition
fh = h^2*fh;
fh(:,1) = fh(:,1) - uh(Iint,1);
fh(:,m) = fh(:,m) - uh(Iint,n+1);
fh(1,:) = fh(1,:) - uh(1,Jint);
fh(m,:) = fh(m,:) - uh(n+1,Jint);

% convert the right hand side from a grid function into a vector
F = reshape(fh,m*m,1);

% Solve the linear system for the grid function uh on the interior, written
% as vector
U = Lh\F;  

% convert from vector into grid function use this for uh at the interior
% grid points (uh is already given at the boundary grid points)
uh(Iint,Jint) = reshape(U,m,m);

% plot results:
figure(1)
clf
uhplot=surf(X,Y,uh);
set(uhplot,'facecolor','yellow')
l=light('position',[-5,5,10]);
lighting gouraud

% skip the rest if the true solution is not known
%return
utrue = u(X,Y);    % utrue is the solution at all the grid points
% print error
err = max(max(abs(uh-utrue)));   
fprintf('Absolute error in max norm = %10.3e\n',err)
fprintf('Relative error in max norm = %6.3f%%\n',100*err/max(max(abs(utrue))))

% plot true solution on fine mesh for comparison
figure(2)
clf
xf = linspace(0,1,101);   % x grid points including boundaries
yf = linspace(0,1,101);   % y grid points including boundaries
[Xf,Yf] = meshgrid(xf,yf);      % 2d arrays of x,y values
Xf = Xf';                     % transpose so that X(i,j),Y(i,j) are
Yf = Yf';                     % coordinates of (i,j) point
uplot=surf(Xf,Yf,u(Xf,Yf));
set(uplot,'edgecolor','none','facecolor','yellow')
l=light('position',[-5,5,10]);
lighting gouraud
