% poisson2.m -- solve the Poisson problem u_{xx} + u_{yy} = f(x,y)
% on [a,b] x [a,b].
%
% The 5-point Laplacian is used at interior grid points.
% This system of equations is then solved using backslash.
%
% From http://www.amath.washington.edu/~rjl/fdmbook/chapter3 (2007)
a = 0;
b = 1;
m = 20;
h = (b-a)/(m+1);
x = linspace(a,b,m+2); % grid points x including boundaries
y = linspace(a,b,m+2); % grid points y 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) point
Iint = 2:m+1; % indices of interior points in x
Jint = 2:m+1; % indices of interior points in y
Xint = X(Iint,Jint); % interior points
Yint = Y(Iint,Jint);
f = @(x,y) 1.25*exp(x+y/2); % f(x,y) function
rhs = f(Xint,Yint); % evaluate f at interior points for right hand side
% rhs is modified below for boundary conditions.
utrue = exp(X+Y/2); % true solution for test problem
% set boundary conditions around edges of usoln array:
usoln = utrue; % use true solution for this test problem
% This sets full array, but only boundary values
% are used below. For a problem where utrue
% is not known, would have to set each edge of
% usoln to the desired Dirichlet boundary values.
% adjust the rhs to include boundary terms:
rhs(:,1) = rhs(:,1) - usoln(Iint,1)/h^2;
rhs(:,m) = rhs(:,m) - usoln(Iint,m+2)/h^2;
rhs(1,:) = rhs(1,:) - usoln(1,Jint)/h^2;
rhs(m,:) = rhs(m,:) - usoln(m+2,Jint)/h^2;
% convert the 2d grid function rhs into a column vector for rhs of system:
F = reshape(rhs,m*m,1);
% form matrix A:
I = speye(m);
e = ones(m,1);
T = spdiags([e -4*e e],[-1 0 1],m,m);
S = spdiags([e e],[-1 1],m,m);
A = (kron(I,T) + kron(S,I)) / h^2;
% Solve the linear system:
uvec = A\F;
% reshape vector solution uvec as a grid function and
% insert this interior solution into usoln for plotting purposes:
% (recall boundary conditions in usoln are already set)
usoln(Iint,Jint) = reshape(uvec,m,m);
% assuming true solution is known and stored in utrue:
err = max(max(abs(usoln-utrue)));
fprintf('Error relative to true solution of PDE = %10.3e \n',err)
% plot results:
clf
hold on
% plot grid:
% plot(X,Y,'g'); plot(X',Y','g')
% plot solution:
contour(X,Y,usoln,30,'k')
axis([a b a b])
daspect([1 1 1])
title('Contour plot of computed solution')
hold off